% non linear lumped system

clear
nn=101;dt=0.02;tc=0.5;

u0=-(1+0.25+(sqrt(2*pi))/15)/tc
figure(1); clf;

for i=1:50
   x(i)=i/50;aa(i)=a(x(i));
end
plot(x,aa, 'linewidth', 1.0)
xlabel('x=T_2^*')
ylabel('a(x)=m_2C_2/m_1C_1')
grid on

 fig=figure(1);
 print(fig,fullfile(pwd, 'Figure6_7.tif'),'-dtiff','-r600')

figure(2); clf;


for n=1:nn
   tn=(n-1)*dt;
   if tn<tc
      u(n)=u0;
   else
      u(n)=0;
   end
   
end
t(1)=0;
yw=[1; 1]  ;y(1)=0.;z(1)=1.;

for n=2:nn
 t(n)=(n-1)*dt; ynew=direct(yw,u(n),dt); 
 y(n)=ynew(1);z(n)=ynew(2);yw=ynew;
  
end

plot(t,y,t,z, 'linewidth', 1.0)
xlabel('t^* = t/\tau = hAt/m_1C_1')
ylabel('T^*')

grid on
legend('T_1','T_2')
text(0.4,0.9,sprintf('t_c = %3.1f',tc) )

% fig=figure(2);
% if( tc == 0.5 )
%     print(fig,fullfile(pwd, 'Figure6_8.tif'),'-dtiff','-r600')
% else
%     print(fig,fullfile(pwd, 'Figure6_9.tif'),'-dtiff','-r600')
% end

function val=a(x)

val=0.25+exp(-225*(x-0.3)^2); 
end

function y=direct(yw,q,dt)
aw=a(yw(2));
A=[1+dt -dt;-dt/aw 1+dt/aw]; 
B=yw;B(1)=B(1)+dt*q;
y=A\B; 
end