bestep(x,h)

The backward Euler method for integrating the system x'=f(x) is

This is an implicit method: given x[i] you need to use a root-finding method, such as Newton's method, to find x[i+1]. Here is a matlab function bestep(x,h) that does exactly this. Put it in a file called bestep.m. The function f defining the system of differential equations should be in the file f.m. f and x should be taken to be column vectors. 

function a=bestep(x,h)

d=size(x,1);

eps=0.00001;
cha=eps*eye(d,d);
M=zeros(d,d);

fx=f(x);
a=x+h*fx;
con=1;

while(con>0)

fa=f(a);

for i=1:d
M(:,i)=(f(a+cha(:,i))-fa)/eps;
end

ch = (eye(d,d)-h*M)\(x+h*fa-a);
con = norm(ch,1);
a=a+ch;

end
Notes: The function uses Newton's method, with derivatives of f computed with the simplest finite difference formula with "small" parameter eps=0.00001. You might want to change this! It also continues to run the Newton iteration until it finds a y such that y=x+hf(y) to machine accuracy. You might want to change this too... just replace the statement "while(con>0)" with "while(con>1e-10)" or something similar.

Back to the course homepage
Back to my main teaching page
Back to my main page