gistfile1.matlab

function bvp( N )

r0 = 0.02; %m
h = 20; %W/m^2*K
k = 14; %W/m*K
TA = 20; %degrees C
T0 = 100; %degrees C
L = 0.1; %m

a = @(x) -1./(L-x/2);
b = @(x) -2*h./(r0*k*(1-x/(2*L)));
c = @(x) 2*h*TA./(r0*k*(1-x/(2*L)));

Dx = L/(N+1/2);


for i=1:N+2
    x(i,1)=(i-1)*Dx;
end

%Create all the other T's

   
sup = (1+Dx/2*a(x(3:N+1)));
main = (b(x(2:N+1))*Dx^2-2);
sub = (1-Dx/2*a(x(2:N)));


A = diag(sub,-1) + diag(main,0) + diag(sup,1);

y = -c(x(2:N+1))*Dx^2;

y(1,1) = y(1,1)-(1-Dx*a(x(1,1))/2)*T0;

A(N,N) = A(N,N)-(1+Dx*a(x(N+2,1))/2)*((h*Dx-2*k)/(h*Dx+2*k));

y(N,1) = y(N,1)-(1+Dx*a(x(N+2,1))/2)*((2*Dx*h*TA)/(h*Dx+2*k));

T = A\y;

Tp1 = (2*Dx*h*TA)/(h*Dx+2*k) - (h*Dx-2*k)/(h*Dx+2*k)*T(N);

TL = (T(N)+Tp1)/2

plot(x(2:N+1),T)


end