I'm trying to solve this but it says ( Attempted to access x(2); index out of bounds because numel(x)=1 ) ..I just want to change the solver from ode32 '' I marked it as a comment'' to modified euler method but it showing this problem

global Pm f H E Y th ngg f=60; R=linedata(:,3); X=linedata(:,4); % added 18/10/2012 %zdd=gendata(:,2)+j*gendata(:,3); ngr=gendata(:,1); %H=gendata(:,4); ngg=length(gendata(:,1)); %% for k=1:ngg zdd(ngr(k))=gendata(k, 2)+j*gendata(k,3); %H(ngr(k))=gendata(k, 4); H(k)=gendata(k,4); % new end %% for k=1:ngg I=conj(S(ngr(k)))/conj(V(ngr(k))); %Ep(ngr(k)) = V(ngr(k))+zdd(ngr(k))*I; %Pm(ngr(k))=real(S(ngr(k))); Ep(k) = V(ngr(k))+zdd(ngr(k))*I; % new Pm(k)=real(S(ngr(k))); % new

end E=abs(Ep); d0=angle(Ep); for k=1:ngg nl(nbr+k) = nbus+k;

nr(nbr+k) = gendata(k, 1);

%R(nbr+k) = gendata(k, 2); %X(nbr+k) = gendata(k, 3);

R(nbr+k) = real(zdd(ngr(k))); X(nbr+k) = imag(zdd(ngr(k)));

Bc(nbr+k) = 0; a(nbr+k) = 1.0; yload(nbus+k)=0; end nbr1=nbr; nbus1=nbus; nbrt=nbr+ngg; nbust=nbus+ngg; linedata=[nl' nr' R X -j*Bc a]; [Ybus, Ybf]=YBUSBF(linedata, yload, nbus1,nbust); fprintf('\nPrefault reduced bus admittance matrix \n') Ybf Y=abs(Ybf); th=angle(Ybf); Pm=zeros(1, ngg); disp([' G(i) E''(i) d0(i) Pm(i)']) for ii = 1:ngg for jj = 1:ngg Pm(ii) = Pm(ii) + E(ii)*E(jj)*Y(ii, jj)*cos(th(ii, jj)-d0(ii)+d0(jj));

end, fprintf(' %g', ngr(ii)), fprintf(' %8.4f',E(ii)), fprintf(' %8.4f', 180/pi*d0(ii)) fprintf(' %8.4f \n',Pm(ii)) end respfl='y'; while respfl =='y' respfl=='Y' nf=input('Enter faulted bus No. -> '); fprintf('\nFaulted reduced bus admittance matrix\n') Ydf=YBUSDF(Ybus, nbus1, nbust, nf) %Fault cleared [Yaf]=YBUSAF(linedata, yload, nbus1,nbust, nbrt); fprintf('\nPostfault reduced bus admittance matrix\n') Yaf resptc='y'; while resptc =='y' resptc=='Y' tc=input('Enter clearing time of fault in sec. tc = '); tf=input('Enter final simulation time in sec. tf = '); clear t x del %%%%%%% from here the changing of solver from ode32 to meuler_simple.%%%%%%% t0 = 0; w0=zeros(1, length(d0)); x0 = [d0, w0]; tol=0.0001; Y=abs(Ydf); th=angle(Ydf);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % tspan=[t0, tc]; % [t1, xf] =ode23('dfpek', tspan, x0);

d= 0.02;N=20; t1=0; h=(tc-t1)/N; for k=1:N t1=t1+h; [t1,x]= meuler_simple(@dfpek, t1,x0,tc,N); x1=x; x1=[x1; x0]; end

x1c =x1(length(x1), :);

....................................... Y=abs(Yaf); th=angle(Yaf); ....................................... % tspan = [tc, tf]; % [t2,xc] =ode23('afpek', tspan, x0c);

l=(tf-tc)/N; t2=tc for n=1:N t2=t2+l; [t2,x]= meuler_simple(@afpek, tc,x1c,tf,N); x2=x; x2=[x2; x1c] end

t =[t1; t2]; x = [x1; x2];