% Lead toxicology problem
% x = concentration in blood
% y = concnetration in soft tissue
% z = concentration in bone
% X = [x;y;z]
% Input = input data = [input,0,0];
% goal solve X'=A*X+Input
% X0 = initial condition =0
function lead
clear all

X0 = zeros(3,1);
Input = [10;0;0];

A = [-0.03610,  0.01240,  0.00004;...
      0.01110, -0.02860,  0.00000 ;... 
      0.00390,  0.00000, -0.00004];

% compute homogeneous solution info
% The diagonal entries of D are the eigenvalues
% the columns of V are the eigenvectors
[V,D] = eig(A);
% To only compute the eigenvalues 
v =eig(A);

% find particular solution info.
Xp = -A\Input;

% now compute coefficients
C = V\(X0-Xp);

% now we can plot the solutions.
t =[0:1:90];

x = C(1)*V(1,1)*exp(D(1,1)*t)+C(2)*V(1,2)*exp(D(2,2)*t)+C(3)*V(1,3)*exp(D(3,3)*t)+Xp(1);
y = C(1)*V(2,1)*exp(D(1,1)*t)+C(2)*V(2,2)*exp(D(2,2)*t)+C(3)*V(2,3)*exp(D(3,3)*t)+Xp(2);
z = C(1)*V(3,1)*exp(D(1,1)*t)+C(2)*V(3,2)*exp(D(2,2)*t)+C(3)*V(3,3)*exp(D(3,3)*t)+Xp(3);

figure(1)
plot(t,x,t,y,t,z)
title('Exact')
xlabel('Time')
ylabel('Concentration')
legend('blood','soft tissue','bone')

[T,Xapp] = ode45(@rhs,[0 90],X0);

figure(2)
plot(T,Xapp(:,1),'b',T,Xapp(:,2),'r',T,Xapp(:,3),'g')
title('Approximate')
xlabel('Time')
ylabel('Concentration')
legend('blood','soft tissue','bone')


return

function dy = rhs(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = -0.03610*y(1)+ 0.01240*y(2)+ 0.00004*y(3) +10;
dy(2) = 0.01110*y(1)+-0.02860*y(2)+ 0.00000*y(3)+0;
dy(3) = 0.00390*y(1)+  0.00000*y(2)+ -0.00004*y(3)+0;

return