% Ritz method - Final Question 1b

clear all;
clc;

syms x T T0 c1 c2 k p h A T_inf L
T = T0+ c1*x+c2*x^2;
dT_dx = diff(T);
F = 0.5*k*dT_dx^2+0.5*(p*h/A)*T^2-(p*h/A)*T_inf*T;
I = int(F,0,L);

e1 = diff(I,c1);
e2 = diff(I,c2);

d = 1e-3; L = 50e-3; T_inf = 30; h = 20; k = 200; T0 = 300;
A = pi*d^2/4; p = pi*d;

e1 = eval(e1);
e2 = eval(e2);

[a b] = solve(e1,e2);

c1 =eval(a);     c2 = eval(b); 

y = [0 10 20 30 40 50]*10^-3;

for i = 1:6
    x = y(i);
    Temp(i) = eval(T);
end
plot(y,Temp,'^:');
legend('Ritz Method');