clc, clear all
disp(['Se rezolva problema variationala cu functionala ce depinde',...
    'de mai multe functii variabile'])
syms x y z Dy D2y Dz D2z % se descriu variabilele simbolice
F=Dy^2-Dz^2+2*y*z-2*y^2; % functia integrand
x1=0; y1=0; z1=0; x2=pi/2; y2=1; z2=1; % conditiile la extremitati
disp('Datele initiale:')
fprintf(['Functia integrand F(x,y,y'',z,z'')=%s\n'],char(F))
disp('Conditiile la extremitatea de stanga:')
fprintf('y(%d)=%d;   z(%d)=%d\n',x1,y1,x1,z1)
disp('Conditiile la extremitatea de dreapta:')
fprintf('y(%d)=%d;   z(%d)=%d\n',x2,y2,x2,z2)
dFdy=diff(F,y); dFdy1=diff(F,Dy); d_dFdy1_dx=diff(dFdy1,x);
d_dFdy1_dy=diff(dFdy1,y); d_dFdy1_dy1=diff(dFdy1,Dy);
d_dFdy1_dz=diff(dFdy1,z); d_dFdy1_dz1=diff(dFdy1,Dz);
dFy1dx=d_dFdy1_dx+d_dFdy1_dy*Dy+...
  d_dFdy1_dy1*D2y+d_dFdy1_dz*Dz+d_dFdy1_dz1*D2z;
dFdz=diff(F,z); dFdz1=diff(F,Dz); d_dFdz1_dx=diff(dFdz1,x);
d_dFdz1_dy=diff(dFdz1,y); d_dFdz1_dy1=diff(dFdz1,Dy);
d_dFdz1_dz=diff(dFdz1,z); d_dFdz1_dz1=diff(dFdz1,Dz);
dFz1dx=d_dFdz1_dx+d_dFdz1_dy*Dy+...
  d_dFdz1_dy1*D2y+d_dFdz1_dz*Dz+d_dFdz1_dz1*D2z;
EulerY=simple(dFdy-dFy1dx); EulerZ=simple(dFdz-dFz1dx);
dEuY=[char(EulerY) '=0']; % EEL in raport cu y
dEuZ=[char(EulerZ) '=0']; % EEL in raport cu z
fprintf('Sistemul EEL:\n%s\n%s\n',dEuY,dEuZ)
Sol=dsolve(dEuY,dEuZ,'x'); % se rezolva sistemul EEL
if length(Sol)~=1
  error('Nu sunt solutii sau sunt mai multe!');
else
  SolY=Sol.y; SolZ=Sol.z;
  disp('Solutia generala a sistemului EEL:')
  fprintf('y(x)=%s\nz(x)=%s\n',char(SolY),char(SolZ))
end
SolLY=subs(SolY,x,x1); % se substituie x1 in y
SolLZ=subs(SolZ,x,x1); % se substituie x1 in z
SolRY=subs(SolY,x,x2); % se substituie x2 in y
SolRZ=subs(SolZ,x,x2); % se substituie x2 in z
EqLY=[char(vpa(SolLY,14)) '=' char(sym(y1))];
EqLZ=[char(vpa(SolLZ,14)) '=' char(sym(z1))];
EqRY=[char(vpa(SolRY,14)) '=' char(sym(y2))];
EqRZ=[char(vpa(SolRZ,14)) '=' char(sym(z2))];
disp('Conditiile la extremitati:')
fprintf('%s\n',EqLY,EqLZ,EqRY,EqRZ)
Con=solve(EqLY,EqLZ,EqRY,EqRZ,'C1,C2,C3,C4');
C1=Con.C1; C2=Con.C2; C3=Con.C3; C4=Con.C4;
Sol3Y=vpa(eval(SolY),14); Sol3Z=vpa(eval(SolZ),14);
disp('Ecuatiile extremalei:')
fprintf('y(x)=%s\nz(x)=%s\n',char(Sol3Y),char(Sol3Z))
Fextr=simple(subs(F,{y,z,Dy,Dz},{Sol3Y,Sol3Z,diff(Sol3Y,x),diff(Sol3Z,x)}));
Jextr=eval(int(Fextr,x,x1,x2))
xpl=linspace(x1,x2); y3=subs(Sol3Y,x,xpl); z3=subs(Sol3Z,x,xpl); 
figure 
plot(xpl,y3,'-r',xpl,z3,'--b') % рисуем график
set(get(gcf,'CurrentAxes'),'FontName','Times New Roman','FontSize',11)
legend('y(x)','z(x)',0)
title('\bf Componentele extremalei admisibile')
xlabel('\itx')
ylabel('\ity\rm(\itx\rm), \itz\rm(\itx\rm)')
grid on, box on
da=daspect; da(1:2)=min(da(1:2)); daspect(da);
figure 
plot3(xpl,y3,z3,'-r') % se deseneaza graficul 3D
set(get(gcf,'CurrentAxes'),'FontName','Times New Roman','FontSize',11)
title('\bfGraficul 3D')
xlabel('\itx')
ylabel('\ity\rm(\itx\rm)')
zlabel('\itz\rm(\itx\rm)')
view(205,30) % se alege punctul de vizualizare
grid on, box on