funcprot(0);

function [yN]=fun(x)
  yN= -(2*real(trace(((expm(-%i*dt*HAMo)'*expm(-%i*dt*(x(1)*ix+x(2)*iy+...
    x(3)*sx+x(4)*sy))- eye(4,4))*A*U*B*A') ))- ...
    dt*l1*((x(1)-x1(j-1)).*(x(1)+x1 (j-1)) +(x(2)-x2(j-1)).*(x(2)+x2(j-1))+...
    (x(3)-x3(j-1)).*(x(3)+x3 (j-1)) +(x(4)-x4 (j-1)).* (x(4)+x4(j-1))));
endfunction
function [d_yN] =d_fun(x)
  d_yN=derivative (fun,x, order = 4);
endfunction
      
function [yN,grd,ind]=f(x,ind)
  yN=fun(x);
  grd=d_fun(x);
endfunction
	   
function [y1,y2,y3,y4]=prop(x1,x2,x3,x4,ix,iy,sx,sy,dt,A,U,B,PROP,l1)
  y1 = x1;
  y2 = x2;
  y3 = x3;
  y4 = x4;
  N=length(x1)+1;
  for j=2:N
    HAMo = rand(4,4);
    B=rand(4,4);
    U=rand(4,4);
    x0= [x1(j-1);x2(j-1);x3(j-1);x4(j-1)];
    printf("BEFORE: %f\n",fun(x0));
    [fopt,xN]=optim(f,x0);
    printf("AFTER: %f\n",fun(xN));
    y1(j-1)	= xN(1);
    y2(j-1)	= xN(2);
    y3(j-1)	= xN(3);
    y4(j-1)	= xN(4);
  end
endfunction

x_pauli=[0,0.5;0.5,0]; 
y_pauli=[0,-0.5*%i;0.5*%i,0]; 
z_pauli=[0.5,0;0,-0.5];
u_pauli=[1,0;0,1];

J = 140;
T = 1/J;  N=10;  NIT=5;
dt=T/N;
rf = 100;
l1 = 1e-4;

//---------------------------------------------------------------------------
ix = x_pauli.*.u_pauli; iy = y_pauli.*.u_pauli; iz = z_pauli.*.u_pauli;
sx = u_pauli.*.x_pauli; sy = u_pauli.*.y_pauli; sz = u_pauli.*.z_pauli;
//---------------------------------------------------------------------------
Hamilt=2*%pi*J*sz*iz;
A=expm(-%i*0.5*dt*Hamilt);

U = rand(4,4);
B = rand(4,4);
PROP=rand(4,4);
//---------------------------------------------------------------------------
u1 = rf*rand(1,N-1); u2 = rf*rand(1,N-1);
u3 = rf*rand(1,N-1); u4 = rf*rand(1,N-1);

[y1,y2,y3,y4]=prop(u1,u2,u3,u4,ix,iy,sx,sy,dt,A,U,B,PROP,l1);

