clear all;
close all;
Delta_t=0.05;
t=0:Delta_t:50;
%
%Preallocating Memory
%
x=zeros(length(t),2);
x(1,1)=1.0;x(1,2)=0.0;
for n=1:length(t)-1
x(n+1,:)=Findxnplus1(Delta_t,x(n,:));
end
%Check with Matlab solution
[tmat,xmat] = ode23(@Example22p4,[0 50],[1 0]);
hold off
plot(t,x(:,1),’k-.’,tmat,xmat(:,1),’k-‘);
hold on
plot(t,x(:,2),’b-.’,tmat,xmat(:,2),’b-‘);
legend(‘x_0(t) Crank-Nicolson’,’x_0(t) MATLAB ode23()’,’x_1(t) Crank-Nicolson’,’x_1(t) MATLAB ode23()’);
xlabel(‘t’,’Fontsize’,18);ylabel(‘x_i(t)’,’Fontsize’,18);
function p=Findxnplus1(Delta_t,xn)
eps=1.0; %Just making sure you go through the loop at least once
p=xn; %Guess value for p
while eps>1.0e-13
[gfunc,jacobian]=Example22p4WithJacobian(p,Delta_t,xn);
eps=max(abs(gfunc));
temp=jacobian\(-gfunc); %this is more efficient than inv(jacobian)*(-gfunc)
p=temp’+p;
end
end
function [gfunc,jac]=Example22p4WithJacobian(p, Delta_t,xn)
gfunc=[p(1)-xn(1)-(Delta_t/2.)*(xn(2)+p(2)); p(2)-xn(2)+(Delta_t/2)*(xn(1)+xn(1).^3+p(1)+p(1).^3)];
jac=[1 -Delta_t/2; (3*p(1).^2+1)*Delta_t/2. 1];
end
function dxdt = Example22p4(t,x)
dxdt = [x(2); -x(1)-x(1).^3];
end