Taylor Method ODE

Question

I am trying to implement the Taylor method for ODEs in MatLab:

My code (so far) looks like this...

function [x,y] = TaylorEDO(f, a, b, n, y0)
% syms t
% x = sym('x(t)'); % x(t)
% f = (t^2)*x+x*(1-x);

h = (b - a)/n;

fprime = diff(f);
f2prime = diff(fprime);

y(0) = y0,
for i=1:n
    T((i-1)*h, y(i-1), n) = double(f((i-1)*h, y(i-1)))+(h/2)*fprime((i-1)*h, y(i-1))
    y(i+1) = w(i) + h*T(t(i), y(i), n); 

I was trying to use symbolic variables, but I don´t know if/when I have to use double.

I also tried this other code, which is from a Matlab function, but I do not understand how f should enter the code and how this df is calculated.

http://www.mathworks.com/matlabcentral/fileexchange/2181-numerical-methods-using-matlab-2e/content/edition2/matlab/chap_9/taylor.m

As error using the function from this link, I got:

>> taylor('f',0,2,0,20)
Error using feval
Undefined function 'df' for input arguments of type 'double'.

Error in TaylorEDO (line 28)
  D = feval('df',tj,yj)

The f I used here was

 syms t
 x = sym('x(t)'); % x(t)
 f = (t^2)*x+x*(1-x);

Answer 1

This is a numerical method, so it needs numerical functions. However, some of them are computed from the derivatives of the function f. For that, you need symbolic differentiation.

Relevant Matlab commands are symfun (create a symbolic function) and matlabFunction (convert a symbolic function to numerical).

The code you have so far doesn't seem salvageable. You need to start somewhere closer to basics, e.g., "Matlab indices begin at 1". So I'll fill the gap (computation of df) in the code you linked to. The comments should explain what is going on.

function [T,Y] = taylor(f,a,b,ya,m)
syms t y
dfs(1) = symfun(f, [t y]);                         % make sure that the function has 2 arguments, even if the user passes in something like 2*y
for k=1:3
    dfs(k+1) = diff(dfs(k),t)+f*diff(dfs(k),y);    % the idea of Taylor method: calculate the higher derivatives of solution from the ODE
end
df = matlabFunction(symfun(dfs,[t y]));            % convert to numerical function; again, make sure it has two variables

h = (b - a)/m;                                     % the rest is unchanged except one line
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m
  tj = T(j);
  yj = Y(j);
  D = df(tj,yj);                                   % syntax change here; feval is unnecessary with the above approach to df
  Y(j+1) = yj + h*(D(1)+h*(D(2)/2+h*(D(3)/6+h*D(4)/24)));
  T(j+1) = a + h*j;
end
end

Example of usage:

syms t y
[T, Y] = taylor(t*y, 0, 1, 2, 100); 
plot(T,Y)