%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)