% Instabilidade numérica

% Solução Exata
alpha = -1;
tau = -1/alpha;
Tmax = 20;
tt = linspace(0,Tmax,10000);
xx = 10*exp(alpha*tt);

% Método de Euler Explícito
figure(1)
clf
x = [];
t = [];
x(1) = 10;
t(1) = 0;
h = 0.25
N = floor(Tmax/h);
for i = 1:N,
  x(i+1) = x(i) + foo(t(i),x(i))*h;
  t(i+1) = t(i) + h;
end

plot(tt,xx,t,x,'*-r')
text(0.75*Tmax,25,sprintf('tau = %3.1f',tau))
xlabel('t')
ylabel('x(t)')
title(sprintf('Estável - h = %3.2f',h))
grid
axis([0,Tmax,-30,30])
shg

figure(2)
clf
x = [];
t = [];
x(1) = 10;
t(1) = 0;
h = 0.5
N = floor(Tmax/h);
for i = 1:N,
  x(i+1) = x(i) + foo(t(i),x(i))*h;
  t(i+1) = t(i) + h;
end

plot(tt,xx,t,x,'*-r')
text(0.75*Tmax,25,sprintf('tau = %3.1f',tau))
xlabel('t')
ylabel('x(t)')
title(sprintf('Estável - h = %3.2f',h))
grid
axis([0,Tmax,-30,30])
shg

figure(3)
clf
x = [];
t = [];
x(1) = 10;
t(1) = 0;
h = 1.0
N = floor(Tmax/h);
for i = 1:N,
  x(i+1) = x(i) + foo(t(i),x(i))*h;
  t(i+1) = t(i) + h;
end

plot(tt,xx,t,x,'*-r')
text(0.75*Tmax,25,sprintf('tau = %3.1f',tau))
xlabel('t')
ylabel('x(t)')
title(sprintf('Estável - h = %3.2f',h))
grid
axis([0,Tmax,-30,30])
shg

figure(4)
clf
x = [];
t = [];
x(1) = 10;
t(1) = 0;
h = 1.5
N = floor(Tmax/h);
for i = 1:N,
  x(i+1) = x(i) + foo(t(i),x(i))*h;
  t(i+1) = t(i) + h;
end

plot(tt,xx,t,x,'*-r')
text(0.75*Tmax,25,sprintf('tau = %3.1f',tau))
xlabel('t')
ylabel('x(t)')
title(sprintf('Estável - h = %3.2f',h))
grid
axis([0,Tmax,-30,30])
shg

figure(5)
clf
x = [];
t = [];
x(1) = 10;
t(1) = 0;
h = 2
N = floor(Tmax/h);
for i = 1:N,
  x(i+1) = x(i) + foo(t(i),x(i))*h;
  t(i+1) = t(i) + h;
end

plot(tt,xx,t,x,'*-r')
text(0.75*Tmax,25,sprintf('tau = %3.1f',tau))
xlabel('t')
ylabel('x(t)')
title(sprintf('Limite do Estável - h = %3.2f',h))
grid
axis([0,Tmax,-30,30])
shg

figure(6)
clf
x = [];
t = [];
x(1) = 10;
t(1) = 0;
h = 2.2
N = floor(Tmax/h);
for i = 1:N,
  x(i+1) = x(i) + foo(t(i),x(i))*h;
  t(i+1) = t(i) + h;
end

plot(tt,xx,t,x,'*-r')
text(0.75*Tmax,25,sprintf('tau = %3.1f',tau))
xlabel('t')
ylabel('x(t)')
title(sprintf('Instável - h = %3.2f',h))
grid
axis([0,Tmax,-30,30])
shg