clear all;
close all;
Delta_t=0.05;
t=0:Delta_t:50;
xRK2=zeros(length(t),3);
xEuler=zeros(length(t),3);
xRK2(1,:)=[0.0 1.0 1.0];
xEuler(1,:)=[0.0 1.0 1.0];
for n=1:length(t)-1
k1=f(t(n),xRK2(n,:))’;
k2=f(t(n)+Delta_t,xRK2(n,:)+Delta_t*k1)’;
xRK2(n+1,:)=xRK2(n,:)+Delta_t*(k1/2.0+k2/2.0);
xEuler(n+1,:)=xEuler(n,:)+Delta_t*f(t(n),xEuler(n,:))’;
end
[tmat,xmat]=ode23(@f,[0 50],[0 1 1]);
hold off
plot(t,xEuler(:,2),’b-‘,’linewidth’,2,’markersize’,10)
hold on
plot(t,xRK2(:,2),’r-‘,’linewidth’,2,’markersize’,10)
plot(tmat,xmat(:,2),’k-‘,’linewidth’,2)
xlabel(‘t’);ylabel(‘x_1(t)’);
legend(‘Explicit Euler’,’RK-2′,’MATLAB ode23()’)
function dxdt=f(t,x)
dxdt=[x(2)*x(3); -x(1)*x(3) ;-0.5*x(1)*x(2)];
end