clear all; close all; % Euler Method t0=0; u0=1; h = input(' Enter the value of the time step ') n=input(' Enter the number of time steps ') t=t0; u=u0; T(1)=t0; sol(1)=u0; for i=1:n k1=f(t,u); u=u+h*k1; t=t+h; T(i+1)=t; sol(i+1)=u; end; % Improved Euler Method t=t0; u=u0; for i=1:n k1=f(t,u); k2=f(t+h,u+h*k1); u=u+0.5*h*(k1+k2); t=t+h; sol_iEm(i+1)=u; end; sol_iEm(1)=u0; % Runge-Kutta Method t=t0; u=u0; for i=1:n k1=f(t,u); k2=f(t+0.5*h,u+0.5*h*k1); k3=f(t+0.5*h,u+0.5*h*k2); k4=f(t+h, u+h*k3); u=u+(h/6)*(k1+2*k2+2*k3+k4); t=t+h; sol_RK(i+1)=u; end; sol_RK(1)=u0; plot(T,sol_exact(T),'r',T,sol,T,sol_iEm,T,sol_RK,'k'); title('Numerical methods for ODEs') legend('exact', 'Euler', 'Modified Euler', 'RK','Location','NorthEast') % The error of these methods for i=1:n+1 E_Euler(i)=abs(sol(i)-sol_exact(T(i))); E_iEm(i)=abs(sol_iEm(i)-sol_exact(T(i))); E_RK(i)=abs(sol_RK(i)-sol_exact(T(i))); end; figure plot(T,E_Euler, T,E_iEm, T, E_RK); title('Error functions'); legend('Euler', 'IEM', 'RK'); %axis ([0 T(n+1) 0 sol(n+1)],'square');