clc; clear all
disp('Se rezolva problema variationala cu functionala ce depinde')
disp('de derivate de ordin superior ale functiei necunoscute')
syms x y Dy D2y D3y D4y % variabilele simbolice
F=D2y^2-16*y^2+x*exp(x);
x1=0; y1=1; Dy1=0; x2=pi/4; y2=0; Dy2=-2;
disp('Datele initiale:')
disp('Functia integrand:')
fprintf('F(x,y,y'',y'''')=%s\n',char(F))
disp('Conditiile la extremitatea de stanga:')
fprintf('y(%d)=%d;   y''(%d)=%d\n',x1,y1,x1,Dy1)
disp('Conditiile la extremitatea de dreapta:')
fprintf('y(%d)=%d;   y''(%d)=%d\n',x2,y2,x2,Dy2)
dFdy=diff(F,y); dFdy1=diff(F,Dy); dFdy2=diff(F,D2y);
d_dFdy1_dx=diff(dFdy1,x); d_dFdy1_dy=diff(dFdy1,y);
d_dFdy1_dy1=diff(dFdy1,Dy); d_dFdy1_dy2=diff(dFdy1,D2y);
dFy1dx=d_dFdy1_dx+d_dFdy1_dy*Dy+d_dFdy1_dy1*D2y+d_dFdy1_dy2*D3y;
d_dFdy2_dx=diff(dFdy2,x); d_dFdy2_dy=diff(dFdy2,y);
d_dFdy2_dy1=diff(dFdy2,Dy); d_dFdy2_dy2=diff(dFdy2,D2y);
dFy2dx=d_dFdy2_dx+d_dFdy2_dy*Dy+d_dFdy2_dy1*D2y+d_dFdy2_dy2*D3y;
d_dFdy2dx_dx=diff(dFy2dx,x); d_dFdy2dx_dy=diff(dFy2dx,y);
d_dFdy2dx_dy1=diff(dFy2dx,Dy); % d((dFy'')/dx)/dy'
d_dFdy2dx_dy2=diff(dFy2dx,D2y); % d((dFy'')/dx)/dy''
d_dFdy2dx_dy3=diff(dFy2dx,D3y); % d((dFy'')/dx)/dy'''
d2Fy2dx2=d_dFdy2dx_dx+d_dFdy2dx_dy*Dy+d_dFdy2dx_dy1*D2y;
d2Fy2dx2=d2Fy2dx2+d_dFdy2dx_dy2*D3y+d_dFdy2dx_dy3*D4y;
Euler=simple(dFdy-dFy1dx+d2Fy2dx2);
deqEuler=[char(Euler) '=0']; % se scrie EEP
fprintf('Ecuatia Euler-Poisson:\n%s\n',deqEuler)
Sol=dsolve(deqEuler,'x'); % se rezolva analitic ecuatia
if length(Sol)~=1 
  error('Nu sunt solutii sau sunt mai multe!');
else
  disp('Solutia generala a EEP:')
  fprintf('y(x)=%s\n',char(Sol))
end
dydx=diff(Sol,x); % derivata
slY=subs(Sol,x,x1); % se substituie x1 in y(x)
slDY=subs(dydx,x,x1); % se substituie x1 in y'(x)
srY=subs(Sol,x,x2); % se substituie x2 in y(x)
srDY=subs(dydx,x,x2); % se substituie x2 in y'(x)
elY=[char(vpa(slY,14)) '=' char(sym(y1))];
elDY=[char(vpa(slDY,14)) '=' char(sym(Dy1))];
erY=[char(vpa(srY,14)) '=' char(sym(y2))];
erDY=[char(vpa(srDY,14)) '=' char(sym(Dy2))];
disp('Conditiile la extremitati:')
fprintf('%s\n',elY,elDY,erY,erDY)
Con=solve(elY,elDY,erY,erDY,'C2,C3,C4,C5');
C2=Con.C2; C3=Con.C3; C4=Con.C4; C5=Con.C5;
Sol4=vpa(eval(Sol),14); % se substituie C2-C5
disp('Ecuatia extremalei:')
fprintf('y(x)=%s\n',char(Sol4))
xpl=linspace(x1,x2); y4=subs(Sol4,x,xpl); 
figure 
plot(xpl,y4,'-r') % se deseneaza graficul
set(get(gcf,'CurrentAxes'),'FontName','Times New Roman','FontSize',11)
title('\bfProblema cu functionala ce contine derivate de ordin superior') 
xlabel('\itx') 
ylabel('\ity\rm(\itx\rm)')
ylim([0 1.1])
da=daspect; da(1:2)=min(da(1:2)); daspect(da);
grid on, box on