Hello,
function euler1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
tic
delta_q=1e-4;
q1=0;q2=0;
x1=0;x2=0;
t=0;delta_t=0;tfinal=20;
A=zeros(0,3);
n=0;
while (t<tfinal)
Dx1=q2;
Dx2=1-3*q1-4*q2;
if (Dx1>0)
delta_x1=delta_q+q1-x1;
else
delta_x1=delta_q-q1+x1;
end
if (Dx2>0)
delta_x2=delta_q+q2-x2;
else
delta_x2=delta_q-q2+x2;
end
delta_t1=delta_x1/abs(Dx1);
delta_t2=delta_x2/abs(Dx2);
if (delta_t1<delta_t2)
delta_t=delta_t1;
t=t+delta_t;
x1=x1+Dx1*delta_t;
x2=x2+Dx2*delta_t;
q1=x1;
else
delta_t=delta_t2;
t=t+delta_t;
x1=x1+Dx1*delta_t;
x2=x2+Dx2*delta_t;
q2=x2;
end
n=n+1;
A(n,1)=t;
A(n,2)=x1;
A(n,3)=x2;
% A(n,4)=q1;
% A(n,5)=q2;
end
figure(1)
close(1)
figure(1)
subplot(1,3,1)
hold on
axis([0,20 -0.05 0.35])
plot(A(:,1),A(:,2),'b')
plot(A(:,1),A(:,3),'r');
grid on
title('HS')
toc
whos
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=linspace(0,20,501);
x0=[0 0];
whos
[t,x]=ode45(@fun,t,x0);
whos
subplot(1,3,2)
hold on
axis([0,20 -0.05 0.35])
plot(t,x(:,1),'k')
plot(t,x(:,2),'k')
grid on
title('edo45')
subplot(1,3,3)
hold on
axis([0,20 -0.05 0.35])
plot(A(:,1),A(:,2),'b')
plot(A(:,1),A(:,3),'r');
plot(t,x(:,1),'k.')
plot(t,x(:,2),'k.')
grid on
title('HS, edo45')
function [dxdt]=fun(~,x)
x1p=x(2);
x2p=1-3*x(1)-4*x(2);
dxdt=[x1p; x2p];
end
end
