Уравнение теплопроводности на конечном стержне

restart;
h:=3; C:=4/5; fo:=3/5; To:=4; u0:=x->To*sin(5*Pi*x/2);
 
                                      3
                                      4
                                      -
                                      5
                                      3
                                      -
                                      5
                                      4
           /5     \
x -> To sin|- Pi x|
           \2     /
 e:=H*cos(k*x)+B*sin(k*x);
 
                           H cos(k x) + B sin(k x)
s2:=eval(subs(x=2,diff(e,x)+h*e));
 
s1:=eval(subs(x=0,diff(e,x)));
 
         -H sin(2 k) k + B cos(2 k) k + 3 H cos(2 k) + 3 B sin(2 k)
                                     B k
B:=solve(s1,B);
                                      0
e:=collect(e,H);
                                 H cos(k x)
f:=collect(s2,H);
                        (-sin(2 k) k + 3 cos(2 k)) H
f:=(k)->(-sin(2*k)*k+3*cos(2*k));
k -> -sin(2 k) k + 3 cos(2 k)
plot(f(k),k=0..20);
 
N:=10;
 
                                     10
for n from 1 to N do 
K[n]:=fsolve(f(k)=0,k=(n-1)*Pi..Pi*n);
H[n]:=1/sqrt(int((3*cos(2*K[n]*x))^2,x=0..2));
end do;
                                0.6747764119
                                0.3601472351
                                 4.983335725
                                0.3299282424
                                 8.032697502
                                0.3307948264
                                 11.12724771
                                0.3317179460
                                 14.24097890
                                0.3322585345
                                 17.36429931
                                0.3325774780
                                 20.49303161
                                0.3327763885
                                 23.62509870
                                0.3329073373
                                 26.75935962
                                0.3329976035
                                 29.89513818
                                0.3330622455
E:=proc(x,n)
global H,K;
return H[n]*(3*cos(2*K[n]*x));
end proc;
proc(x, n)  ...  end;
plot([E(x,1),E(x,2),E(x,3)],x=0..2);
 
plot(u0(x), x=0..2);
 
for n from 1 to N do #ðàñêëàäûâàåì íàøå íà÷àëüíîå óñëîâèìå ïî áàçèñó
G[n]:=int(E(x,n)*u0(x),x=0..2);
end do;
                                0.05460857023
                                -1.212419451
                                -0.2786021127
                                -0.1344403407
                               -0.08000677208
                               -0.05319198572
                               -0.03795072495
                               -0.02844712498
                               -0.02211879790
                               -0.01769184138
T:=proc(x,t)
global N,E,G;local n;
return sum(G[n]*E(x,n)*exp(-K[n]^2*t),n=1..N);
end proc;
proc(x, t)  ...  end;
plot3d(T(x,t),x=0..1,t=0.1..1);