%Matlab script to calculate temperature distribution in
%a fin using the Shooting method. Here we convert a bounary
%value problem into an initial value problem and use a
%4th order Runge-Kutta solver.
%Note that k1,k2,k3 and k4 are different slope
%estimates.
%Calls the function Derivs which calculates the derivative
%of the function.
%
%In this example there are two first order ODEs that must be solved
%with the initial conditions y(1) and y(2)
clear;
clc;

x=0;            %initial x- condition

y(1)=4;            %initial y1-condition - T(0)
y(2)=0;         %intital y2-condition - Guess dT/dx(0) : -10

xnew(1)=x;
ynew(1,1)=y(1);
ynew(1,2)=y(2);

m=length(y);     %number of equations

h = 0.25;        %step size
xend = 2;       %upper x-range of integration
n=(xend-x)/h;   %number of integration steps

for i=2:n+1
    
%calculate k1 and 1st intermdiate y value
k1=derivs_multi(x,y)          % k1 - slope at (xi,yi)
ym=y + k1*h/2;                 % ym - 1st intermediate y value

%calculate k2 and 2nd intermediate y-value based on k2
k2=derivs_multi(x+h/2,ym)    % k2 - slope at (xi+h,y+k1*h/2)
ym=y + k2*h/2;                % ym - 2nd intermediate y value

%calculate k3 and 3rd intermediate y-value based on k2
k3=derivs_multi(x+h/2,ym)     % k3 - slope at (xi+h/2,y+k2*h/2)
ye=y + k3*h;                   % ye - last y value

%calculate k4 
k4=derivs_multi(x+h,ye)      % k4 - slope at (x+h,y+k3*h)

%calculate the weighted average slope
slope=(k1 + 2*(k2 + k3) + k4)/6
%calculate y value at new x-location
ynew(i,:)=y+slope*h;
y=y+slope*h;
xnew(i)=x+h;
x=x+h;
end

%send the integrated values to the screen
xnew
ynew
%Get Exact Solution and plot with RK4 solution
% yexact = 3.881524*exp(sqrt(0.1)*xnew) + 196.1185*exp(-sqrt(0.1)*xnew);
yexact = exp(-xnew/4).*(4*cos(2.6339*xnew) + 0.37966*sin(2.6339*xnew)) 
plot(xnew,ynew(:,1),'*',xnew,yexact);
key1string = ['RK4 solution y2_{guess} = ',num2str(ynew(1,2))];
legend(key1string,'Exact Solution');
xlabel('x');ylabel('y');


%This Section is for linear ODEs only! and should be modified based on 
%running the code twice
%Need to linearly interpolate the 2 solutions to derive the appropriate
%initial condition
tact = 100;         %The actual temperature at right hand bounary
y2guess1 = -10;
y2guess2 = -100;
tguess1 = 1993.4;
tguess2 = -1361.7;
y2guess_final = y2guess1 + ((y2guess2 - y2guess1)/(tguess2 - tguess1))*(tact-tguess1)