ODE45

y0=1;
xspan=[0,1];
f=@(x,y)3*x+y/2;
[x,y]=ode45(f,xspan,y0);
plot(x,y,'*')
        

Copy Code

dsolve

y=dsolve('Dy=3*x+y/2','y(0)=1','x')
x=0:0.1:1;
z=eval(y);
plot(x,z,'*')
        

Copy Code

polyfit

x=[ 10 12 16 20 25 30];
y=[29 33 41 53 65 70];
p=polyfit(x,y,1);
x1=linspace(x(1),x(end));
y1=polyval(p,x1);
plot(x,y,'*',x1,y1)
legend('data points','linear curve')
        

Copy Code

EULARS MODIFIED METHOD

function [x,y] = modEulerODE(x0, xn, y0, h, f)
x=x0:h:xn;
n=length(x)-1;
y=zeros(1,n+1);
y(1) = y0;
for j=1:n
 yp = y(j) + h*f(x(j), y(j));
 y(j+1) = y(j) + 0.5*h*( f(x(j), y(j)) + f(x(j+1), yp));
end
fprintf('Approximate solution of the ODE is\n');
disp([x;y]);
ye=eval(dsolve('Dy=y-x^2+1','y(0)=0.5','x'));% Change according to the problem
plot(x,y,'*',x,ye,'r');
end
        

Copy Code

MILNE'S METHOD

function [x,y]=MilneODE(x0,xn, y0,h,f )
x=x0:h:xn;
n=length(x)-1;
y = zeros(1,n+1);
y(1)=y0;
for i=1:3
 K1=h*f(x(i),y(i));
 K2=h*f(x(i)+0.5*h,y(i)+0.5*K1);
 K3=h*f(x(i)+0.5*h,y(i)+0.5*K2);
 K4=h*f(x(i)+h,y(i)+K3);
 y(i+1)=y(i)+(1/6)*(K1+2*K2+2*K3+K4);
end
for i=4:n
 y(i+1)=y(i-3)+(4*h/3)*(2*f(x(i),y(i)) - f(x(i-1),y(i-1))+2*f(x(i-2),y(i-2)));
 y(i+1)=y(i-1)+(h/3)*(f(x(i+1),y(i+1)) + 4*f(x(i),y(i)) + f(x(i-1),y(i-1)));

end
fprintf('the appx solution of ode is/n')
disp([x;y])
ye=eval(dsolve('Dy=y-x^2+1','y(0)=0.5','x'));
plot(x,y,'*',x,ye)
end
        

Copy Code

Adams Bash forth Method

function [x,y]=Adams_method(x0,y0,xn,n)
h=(xn-x0)/n;
f=@(x,y) x^2*(1+y);
x=x0:h:xn;
y=zeros(1,n);
maxit=10;
s(:,:)=zeros(n,maxit);
y(1)=y0;
for i=1:3
    K1=h*f(x(i),y(i));
    K2=h*f(x(i)+0.5*h,y(i)+0.5*K1);
    K3=h*f(x(i)+0.5*h,y(i)+0.5*K2);
    K4=h*f(x(i)+h,y(i)+K3);
    y(i+1)=y(i)+(1/6)*(K1+2*K2+2*K3+K4);
end
for i=4:n
    y(i+1)=y(i)+(h/24)*(55*f(x(i),y(i))-...
        59*f(x(i-1),y(i-1))+37*f(x(i-2),y(i-2))-...
        9*f(x(i-3),y(i-3)));
    s(i+1,1)=y(i+1);
    for j=1:maxit
        s(i+1,j+1)=y(i)+(h/24)*(f(x(i-2),y(i-2))-...
        5*f(x(i-1),y(i-1))+19*f(x(i),y(i))...
        +9*f(x(i+1),s(i+1,j)));
        if (s(i+1,j+1)-s(i+1,j))<10^(-6)
            y(i+1)=s(i+1,j+1);
               break
        end
    end
end
plot(x,y,'*')
hold on
u=eval(dsolve('Dy=0.09*y-20000','y(0)=200000','x'));
plot(x,u,'r')
end
        

Copy Code

EXPLICIT HEAT METHOD

function oneD_heat(t0, tn,x0,xn,h,k,c )
t=t0:k:tn;
x=x0:h:xn;
m=length(x);
n=length(t);
a=c*k/h^2;
f=@(x) sin(pi*x);
u=zeros(m,n);
u(:,1)=f(x);
if a> 0.5
 fprintf('The method fails')
 return
end
for j=1:n-1
 for i=2:m-1
 u(i, j+1)=a*u(i-1,j)+(1-2*a)*u(i,j)+a*u(i+1,j);
 end
end
for j=1:n
 plot(x,u(:,j))
 hold on
end
figure
surf(t,x,u)
xlabel('t')
ylabel('x')
zlabel('u')
end
        

Copy Code

RUNGE KUTTA METHOD

function [x,y]=RK4ODE(x0,xn, y0,h,f )
x = x0:h:xn;
n = length(x)-1;
y = zeros(1,n+1);
y(1) = y0;
for i=1:n
 K1 = h*f(x(i),y(i))
 K2 = h*f(x(i)+0.5*h,y(i)+0.5*K1)
 K3 = h*f(x(i)+0.5*h,y(i)+0.5*K2)
 K4 = h*f(x(i)+h,y(i)+K3)
 y(i+1) = y(i)+(1/6)*(K1+2*K2+2*K3+K4)
end
fprintf('the appx solution of ode is/n')
disp([x;y])
ye=eval(dsolve('Dy=y-x^2+1','y(0)=0.5','x'));
plot(x,y,'*',x,ye)
end
        

Copy Code