ppp

Euler’s Method MATLAB Program. 

function[]=m_euler 

h=0.1; 

f=@(x,y) x.^2+y; 

x=0:h:1; 

y(1)=1; 

n=length(x); 

yp(1)=y; 

maxit=10; 

for i=2:n 

    yp(i)=y(i-1)+h*f(x(i-1),y(i-1)); 

    s(i,1)=yp(i); 

    for j=2:maxit 

        s(i,j) = y(i-1)+(h/2)*(f(x(i-1),y(i-1))+f(x(i),s(i,j-1))); 

        if abs(s(i,j)-s(i,j-1))<10^(-6)

 

 

            y(i)=s(i,j); 

            break 

        end 

    end 

end 

fprintf('Appr. sol. y(%f)=%f',x(end),y(end)) 

plot(x,yp,'m--') 

hold on 

plot(x,y,'r--') 

hold on 

u=dsolve('Dy=x^2+y','y(0)=1','x'); 

u=eval(u); 

plot(x,u,'b--') 

legend('Appr. sol. by Euler','Appr. sol. by Modified Euler','Exact sol.') 

end

dsolve and ode45

y=dsolve('Dy=3*x+y/2','y(0)=1','x')

x=0:0.1:1;

z=eval(y)

plot(x,z,'X')

hold on

y0=1;

xspan=[0,1];

f=@(x,y) 3*x+y/2;

[x,y]=ode45(f,xspan,y0)

plot(x,y,'r','LineWidth',2)

oneD_heat3 

function[]= oneD_heat3 

t0=0; 

tn=0.5; 

x0=0; 

xn=1; 

h=0.1; 

k=0.005; 

c=1; 

t=t0:k:tn; 

x=x0:h:xn; 

m=length(x); 

n=length(t); 

a=c*k/h^2; 

f=@(x) x.*(1-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 

surf(t,x,u)

xlabel('t') 

ylabel('x') 

zlabel('u') 

end

milne_RK4