reach to phase portrait in a 4*4 matrix and plotting all of them in one figure (if possible), otherwise plotting them 2 by 2 versus to each other.

Hi guys. I wanna draw phase portarait of a matrix 4*4 that is solved for each of its 2 states based on each other. when I write this code, it gives me this error:

Error in

@(t,V)[(3.7*10^-15*V(1)+7.6*10^-6*V(2));(-4*10^-16*V(1)-6.5*10^-5*V(2));(2.37*10^-16*V(1)-0.142*V(3));(-8.15*10^-14*V(1)-0.5535*V(2)+18.21*V(3)+0.0833*V(4))]

Error in test111 (line 70)

Yprime = f1(t,[x(j); y(j)]);

And also why when I change the initial conditions to ones which I comment, it does not give me the same solution after 2500 iteration? Actually it is the operation point.(x1(1)=57.5,x2(1)=8.5 x3(1)=0.051, x4(1)=4.8) and I expect that I reach around this states whatever initial condition that I start.

Thanks

clc;
clear all;
close all;
x1(1)=57.5;
x2(1)=8.5;
x3(1)=0.051;
x4(1)=4.8;
% x1(1)=80;
% x2(1)=10;
% x3(1)=1;
% x4(1)=8;
deltat=10^-3;
A1=[3.4*10^-15 7.6*10^-6 0 0;-4.08*10^-16 -6.55*10^-5 0 0;2.37*10^-16 0.00059 -0.14 0;-8.15*10^-14 -0.05 18.216 0.083];
A2=[9.04*10^-7 -1.48*10^-6 0 0;-1.97*10^-7 -1*10^-5 0 0;9.24*10^-8 0.000458 -0.12 0;-8.55*10^-5 -0.104 65.217 -0.083];
B1=[0.0015 -0.0015 -6.9*10^-12;-5.9*10^-5 -9*10^2 5.9*10^-11;3.46*10^-5 5.2*10^-5 3.24*10^-11;-0.0167 0.0239 -1.0733*10^-9];
B2=[0.0014 -0.0013 -1.9*10^-11;-2.9*10^-5 -0.00011 6.1*10^-11;1.36*10^-5 5.33*10^-5 2.15*10^-11;-0.01 0.042 -4.47*10^-8];
u=[32;32;80];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P=[161469.069263159,8918967.27830892,-576.987918162128,1.43011607072035;8918967.27830892,494939491.276476,-3069582.33611650,-8755.93015894571;-576.987918162128,-3069582.33611650,708027606.278748,206347.846527781;1.43011607072035,-8755.93015894571,206347.846527781,4552.55710715823];
for i=1:2500
% x1(i+1)=(3.7*10^-15*x1(i)+7.6*10^-6*x2(i)+0.554*10^-9)*deltat+x1(i);
% x2(i+1)=(-4*10^-16*x1(i)-6.5*10^-5*x2(i)-2.8)*deltat+x2(i);
% x3(i+1)=(2.37*10^-16*x1(i)-0.142*x3(i)-0.16*10^-2+0.275*10^-2)*deltat+x3(i);
% x4(i+1)=(-8.15*10^-14*x1(i)-0.5535*x2(i)+18.21*x3(i)+0.0833*x4(i)+0.2)*deltat+x4(i);
% V1(i)=0.5*(x1(i).^2+x2(i).^2+x3(i).^2+x4(i).^2);
x1(i+1)=(9.04*10^-7*x1(i)-1.48*10^-6*x2(i)+0.0032)*deltat+x1(i);
x2(i+1)=(-1.97*10^-7*x1(i)-1.008*10^-5*x2(i)-0.0044)*deltat+x2(i);
x3(i+1)=(9.2*10^-8*x1(i)+0.000458*x2(i)-0.127*x3(i)+0.0021)*deltat+x3(i);
x4(i+1)=(-8.55*10^-5*x1(i)-0.1*x2(i)+65.217*x3(i)-0.0833*x4(i)+1.024)*deltat+x4(i);
V2(i)=[x1(i) x2(i) x3(i) x4(i)]*P*[x1(i) x2(i) x3(i) x4(i)]';
%V2(i)=x1(i).^2+300000000*x2(i).^2+200000000*x3(i).^2+0.00000001*x4(i).^2;
end 

for i=1:2500
    figure(1)
plot(x1(1,:),x2(1,:))
hold on;
    plot(x1(1,1),x2(1,1),'bo') % starting point
    hold on;
     plot(x1(1,end),x2(1,end),'ks') % ending point 
     hold on;
     figure(2);
     plot(i,x1(1,i),'*-')
     hold on;
%  figure(3)
%  plot(i,V1(1,i),'*-')
%  hold on
figure(4) 
 plot(i,V2(1,i),'*-')
hold on;
end
hold off;
%---------------------- Vector field & phase portrait
f1=@(t,V)[(3.7*10^-15*V(1)+7.6*10^-6*V(2));(-4*10^-16*V(1)-6.5*10^-5*V(2));(2.37*10^-16*V(1)-0.142*V(3));(-8.15*10^-14*V(1)-0.5535*V(2)+18.21*V(3)+0.0833*V(4))];
f2=@(t,V)[(9.04*10^-7*V(1)-1.48*10^-6*V(2));(-1.97*10^-7*V(1)-1.008*10^-5*V(2));(9.2*10^-8*V(1)+0.000458*V(2)-0.127*V(3));(-8.55*10^-5*V(1)-0.1*V(2)+65.217*V(3)-0.0833*V(4))];
v1 = linspace(-10,10,50);
v2 = linspace(-10,10,50);
v3 = linspace(-10,10,50);
v4 = linspace(-10,10,50);
[x,y] = meshgrid(v1,v2);
size(x)
size(y)
u = zeros(size(x));
v = zeros(size(y));
s = zeros(size(x));
z = zeros(size(y));
t=0; % we want the derivatives at each point at t=0, i.e. the starting time
for j = 1:numel(x)
    Yprime = f1(t,[x(j); y(j)]);
    Yprim = f2(t,[x(j); y(j)]);
    u(j) = Yprime(1);
    v(j) = Yprime(2);
    s(j)= Yprim(1);
   z(j) = Yprim(2);
end
 figure(5)
quiver(x,y,u,v,'r');

figure(6)
 quiver(x,y,s,z,'b');
xlabel('x1')
ylabel('x2')
legend('vector field')
axis tight equal;

 for y10=[0 1]
for y20=[0 1]
   [ts,ys] = ode45(f1,[0,200],[y10;y20]);
   [ts1,ys1] = ode45(f2,[200,350],[y10;y20]);
   [ts2,ys2] = ode45(f1,[350,500],[y10;y20]);
   [ts3,ys3] = ode45(f2,[500,700],[y10;y20]);
   [ts4,ys4] = ode45(f1,[700,900],[y10;y20]);
   [ts5,ys5] = ode45(f1,[900,1000],[y10;y20]);
   figure(7);
   plot(ys(:,1),ys(:,2))
    hold on;
       plot(ys(1,1),ys(1,2),'bo') % starting point
    hold on;
     plot(ys(end,1),ys(end,2),'ks') % ending point
     hold on;
      plot(ys1(:,1),ys1(:,2))
        hold on;
     ylim([-10 100])
    plot(ys1(1,1),ys1(1,2),'bo') % starting point
     hold on;
     plot(ys1(end,1),ys1(end,2),'ks') % ending point  
       hold on;
    plot(ys2(:,1),ys2(:,2))
      hold on;
    plot(ys3(:,1),ys3(:,2))
      hold on;
    plot(ys4(:,1),ys4(:,2))
 hold on;
    plot(ys5(:,1),ys5(:,2))
    xlabel('x1')
    ylabel('x2')
    legend('Phase Portrate x1,x2')
    end
end