% Ritz method - Final Question 1b

clear all;
clc;

d = 1e-3; L = 50e-3; T_inf = 30; h = 20; k = 200; T0 = 300;

A = pi*d^2/4; p = pi*d;

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

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

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

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

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

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