Plot error: index in position 1 is invalid

I want to plot Psi(0,t) and save it as a figure (line 285) in the long code below but got error: Index in position 1 is invalid. Array indices must be positive integers or logical values. How to fix this please? The code is also attached.

% colP3D
clear; clf;
% penetration flag 0 = PKN, 1 = P3D, dimensionless: toughness, leak-off;
ipen = 1;                                      K = 0;      C = 1;
% # of Colloc pts, Geometric time factor,  Max Newton iters,    residual Tol
N = 30;           rfac = 1.2;                Maxit = 30;           tol = 1e-9;
% some constants
gamma_0 = 1.0006328466775270;
gamma_mt = 0;
t0 = 1e-6;
tM = 50;
tc = tM;
if C > 0
    gamma_mt = 2/pi/C;
    Cp = 10/3;
    tc = (gamma_mt/gamma_0)^Cp;
    t0 = 1e-6/C^Cp;
    tM = 1e5/C^Cp;
end
%
lambda_u    = (8+sqrt(K^4+64))/K^2; % runaway height growth boundary

% set up mesh and time steps
t = t0;
dt = t0/10;
dx = 1/N;
x = 0:dx:1;
xf = 0:dx/20:1;
xft = 1-dx:dx/50:1;
in = 1:N;
xi = x(in);
xh = 0.5*(x(1:N)+x(2:N+1));
[xph,is] =  sort([x xh(1:N-1)]);

% precalculate PHI(OMEGA)
[PHI,LAMBDA,OMEGA] = PHIOMEGA(N,K,lambda_u,ipen);

% initialize the solution to the storage PKN values
[Omega0,P0,lambda0,Psi0,gamma0,gammadot] = init_PKN(x,t0); %#ok<ASGLU> % collocation pts
[Omega0h,P0h,lambda0h,Psi0h,gamma0,gammadot] = init_PKN(xh,t0); % 1/2 pts

% now initialize the collocation soln at channel points
y0 = [Omega0(in);
    Psi0(in)];
y = y0;
y0h = [Omega0h(in);
    Psi0h(in)];
lambda = ones(N+1,1);
itime = 0;
tauh = 0:t/10:t;
gammah = gamma_0*tauh.^(4/5);
gamma = gamma_0*(t+dt).^(4/5);
xi_l = 0;
fail = false;
while t < tM && lambda(1) < min(5,lambda_u)
    t = t + dt;
    if t >=  tM; break; end
    R = [];
    tauh = [tauh t]; %#ok<AGROW>
    gammah = [gammah gamma]; %#ok<AGROW>
    for iter = 1:Maxit
        b = getRHS(xi,y,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
        J =  getJacN(xi,y,gamma,K,t,tc,y0,y0h,gamma0,dt,b,OMEGA,PHI,C,tauh,gammah,ipen);
        dS = J\b;
        y = y-reshape(dS(1:2*N),2,N);
        gamma = gamma-dS(2*N+1);
        gammah(end) = gamma;
        Omega = [y(1,:) 0];
        lambda = LamEval(OMEGA,LAMBDA,Omega);
        if lambda(1) > lambda_u
            disp(' lambda_u exceeded ');
            break;
        end
        Psi = [y(2,:) 0];
        if iter > 2 && norm(b,2) < tol; break; end
    end
    if norm(b,2) >= tol
        fprintf('No convergence at t = %g out of %g\n', t, tM);
        fail = true;
        break;
    end
    [yh,Vol,Leak,Atau,alp] = getyh(xi,y,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
    lambda = LamEval(OMEGA,LAMBDA,Omega);
    itime = itime+1;
    TS(itime) = t; %#ok<SAGROW>
    GAMS(itime) = gamma; %#ok<SAGROW>
    ip = find(lambda > 1,1,'last');
    if ip > 1
        xi_l = min(x(ip-1)+(1-lambda(ip-1))/(lambda(ip)-lambda(ip-1))*(x(ip)-x(ip-1)),1);
    end
    XIL(itime)             = xi_l; %#ok<SAGROW>
    plotsol(N,x,Omega,xf,t,K,C,Psi,ipen,xft,Atau,alp,lambda,lambda_u,TS,GAMS,gamma_0,gamma_mt,XIL)
    % copy vectors over for the next time step
    y0 = y;
    y0h = yh;
    dgamma = gamma-gamma0;
    gamma0 = gamma;
    gamma = gamma0+dgamma;
    dt =  rfac*dt;
end
if fail
    disp('Stopped early due to failure');
else
    disp('All done!');
end

function [PHI,LAMBDA,OMEGA] = PHIOMEGA(N,K,lambda_u,ipen)
    % set up form factor PHI a priori  at sample points ready for spline interp
    if ipen == 0
        OMEGA = (1:N);
        LAMBDA = ones(1,N);
        PHI =  ones(1,N);
    else
        Omegafun = @(x,k) (k*x.^(3/2)+2*sqrt(x.^2-1))/pi;
        if K == 0
            lambdaM = 10;
        else
            lambdaM = lambda_u;
        end
        Ns = 1500;
        ns = 1:Ns;
        a = 1.11;
        ep = (exp(lambdaM/Ns/a)-1);
        LAMBDA = [1:0.001:1.1 a*(1+ep).^ns];
        OMEGA = Omegafun(LAMBDA,K);
        PHI = Phi(LAMBDA,K);
    end
end

function y = Phi(lambda,K)
    y = ones(1,length(lambda));
    i2 = find(lambda > 1);
    lam  = lambda(i2);
    y(i2) = pi^2*(4*sqrt(lam)-K*sqrt(lam.^2-1)).*IntOmega3(lam,K)./(lam.^2.*(4*sqrt(lam)+3*K*sqrt(lam.^2-1)))./(((K*lam.^(3/2)+2*sqrt(lam.^2-1))/pi).^3)/12;
end

function y = IntOmega3(lambda,K)
    y = zeros(1,length(lambda));
    for i = 1:length(lambda)
        lam = lambda(i);
        y(i) = 2*quadgk(@(x)OmegaF3(x,lam,K),0,0.5*lam,'abstol',1e-12,'reltol',1e-8);
    end
end

function Omega3 = OmegaF3(x,lambda,K)
    Omega3 = ((4/pi)*(((K/pi/sqrt(lambda)-(2/pi)*asin(1/lambda))*sqrt(lambda^2-4*x.^2))+(2/pi)*(sqrt(lambda^2-4*x.^2)*asin(1/lambda)+phi12(x,lambda)))).^3;
end

function y = phi12(x,lambda)
    i0 = abs(x) == 0.5;
    i1 = find(abs(x) < 0.5);
    i2 = find(0.5 < abs(x) & abs(x) <= 0.5*lambda);
    y(i0) = log(lambda);
    y(i1) = -2*x(i1).*atanh((2*x(i1)*sqrt(lambda^2-1))./sqrt(lambda^2-4*x(i1).^2))+atanh((sqrt(lambda^2-1))./sqrt(lambda^2-4*x(i1).^2));
    y(i2) = -2*x(i2).*atanh(sqrt(lambda^2-4*x(i2).^2)./(2*x(i2)*sqrt(lambda^2-1)))+atanh(sqrt(lambda^2-4*x(i2).^2)./sqrt(lambda^2-1));
end

function [Omega0,P0,lambda0,Psi0,gamma0,gammadot] = init_PKN(x,t0)
    % initialize the solution at collocation points to the impermeable PKN
    gamma_0 = 1.0006328466775270;
    Omega_0 = ((12/5)*gamma_0^2)^(1/3);
    N = length(x);
    gamma0 = gamma_0*t0^(4/5);
    gammadot = (4/5)*gamma_0*t0^(-1/5);
    Omega0 = Omega_0*t0^(1/5)*(1-x).^(1/3).*(1-(1-x)/96 +23*(1-x).^2/64512);
    P0 = Omega0;
    lambda0 = ones(N,1);
    Psi0 = t0^(1/9)*(1-x).^(1/3);
end

function b = getRHS(x,y,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen)
    [n,N] =  size(y); %#ok<ASGLU>
    F  = FS(x,y,gamma,K,y0,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
    h = diff(x);
    H = spdiags(h',0,N-1,N-1);
    ii = 1:N-1;
    ip1 = 2:N;
    xi = x(ii);
    yi = y(:,ii);
    Fi = F(:,ii);
    xip1 = x(ip1);
    yip1 = y(:,ip1);
    Fip1 = F(:,ip1);
    xih = 0.5*(xi+xip1);
    yih = 0.5*(yi+yip1)-(Fip1-Fi)*H/8;
    Fih = FS(xih,yih,gamma,K,y0h(:,ii),gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
    Phim = yip1-yi-(Fip1+4*Fih+Fi)*H/6;
    dgamma = (gamma-gamma0)/dt; %
    if t < tc
        alp = 1/3;
        Atau = (3*dgamma*gamma)^(1/3);
    else
        alp =  3/8;
        Atau = 2*(C/3)^(1/4)*gamma^(3/8)*dgamma^(1/8);
    end
    Voltip = Atau*h(1)^(1+alp)/(1+alp);
    tipVal = Atau*h(1)^alp;
    Vol = (yip1(1,:)+4*yih(1,:)+yi(1,:))*h'/6+Voltip;
    Leaktip = (2/3)*(gamma*dt/(gamma-gamma0))^(1/2)*h(1)^(3/2);
    Leak = (sqrt(t-spline(gammah,tauh,gamma*xip1))+4*sqrt(t-spline(gammah,tauh,gamma*xih))+sqrt(t-spline(gammah,tauh,gamma*xi)))*h'/6+Leaktip;
    Phig = gamma*(Vol+2*C*Leak) - t;
    Phib = GS(y(:,1), y(:,N), tipVal,t);
    b = [Phib(:);
        Phim(:);
        Phig];
end

function F = FS(x,y,gamma,K,y0,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen) %#ok<INUSL>
    Ups = Upsilon(OMEGA,PHI,y(1,:),ipen);
    tau = tauh(end);
    F(1,:)  = -gamma*y(2,:)./y(1,:).^3./Ups;
    tau0 = spline(gammah,tauh,gamma*x);
    F(2,:)  = -gamma*(y(1,:)-y0(1,:))/dt-(gamma-gamma0)*gamma*x.*y(2,:)./y(1,:).^3./Ups/dt-C*gamma./sqrt(tau-tau0);
end

function Ups = Upsilon(OMEGA,PHI,Omega,ipen)
    if ipen == 0
        Ups = ones(1,length(Omega));
    else
        i1 = find(Omega <= OMEGA(1));
        i2 = find(Omega > OMEGA(1));
        Ups(i1) = ones(1,length(i1));
        Ups(i2) = spline(OMEGA,PHI,Omega(i2));
    end %
end

function r = GS(ya, yb, tipVal,t)
    alpha=1/9;
    r(1) = ya(2)-t.^alpha;
    r(2) = yb(1) - tipVal;
end

function J = getJacN(x,y,gamma,K,t,tc,y0,y0h,gamma0,dt,b,OMEGA,PHI,C,tauh,gammah,ipen)
    ep = sqrt(eps);
    [n,N] = size(y);
    I = eye(n*N);
    NC = 1:(n*N+1);
    J = zeros(n*N+1, n*N);
    for j = 1:n*N
        yt = reshape(y(:)+ep*I(:,j),n,N);
        J(NC,j) = (getRHS(x,yt,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen)-b)/ep;
    end
    J(NC,n*N+1) = (getRHS(x,y,gamma+ep,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,[gammah(1:end-1) gammah(end)+ep],ipen)-b)/ep;
end

function lambda = LamEval(OMEGA,LAMBDA,Omega)
    i1 = find(Omega <= OMEGA(1));
    i2 = find(Omega > OMEGA(1));
    lambda(i1) = ones(1,length(i1));
    lambda(i2) = spline(OMEGA,LAMBDA,Omega(i2));
end

function [yih,Vol,Leak,Atau,alp] = getyh(x,y,gamma,K,t,tc,y0,y0h,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen) %#ok<INUSL>
    [n,N] = size(y); %#ok<ASGLU>
    F = FS(x,y,gamma,K,y0,gamma0,dt,OMEGA,PHI,C,tauh,gammah,ipen);
    h = diff(x);
    H = spdiags(h',0,N-1,N-1);
    ii = 1:N-1;
    ip1 = 2:N;
    xi = x(ii);
    yi = y(:,ii);
    Fi = F(:,ii);
    xip1 = x(ip1);
    yip1 = y(:,ip1);
    Fip1 = F(:,ip1);
    xih  = 0.5*(xi+xip1);
    yih = 0.5*(yi+yip1)-(Fip1-Fi)*H/8;
    dgamma = (gamma-gamma0)/dt;
    if t < tc
        alp = 1/3;
        Atau = (3*dgamma*gamma)^(1/3);
    else
        alp =  3/8;
        Atau = 2*(C/3)^(1/4)*gamma^(3/8)*dgamma^(1/8);
    end
    Voltip = Atau*h(1)^(1+alp)/(1+alp);
    tipVal =  Atau*h(1)^alp; %#ok<NASGU>
    Vol = (yip1(1,:)+4*yih(1,:)+yi(1,:))*h'/6+Voltip;
    Leaktip = (2/3)*(gamma*dt/(gamma-gamma0))^(1/2)*h(1)^(3/2);
    Leak = (sqrt(t-spline(gammah,tauh,gamma*xip1))+4*sqrt(t-spline(gammah,tauh,gamma*xih))+sqrt(t-spline(gammah,tauh,gamma*xi)))*h'/6+Leaktip;
end

function plotsol(N,x,Omega,xf,t,K,C,Psi,ipen,xft,Atau,alp,lambda,lambda_u,TS,GAMS,gamma_0,gamma_mt,XIL)
figure(4);
plot(t,Psi(0,t),'bo--','linewidth',2,'markerfacecolor','b')
xlabel('\tau');
ylabel('\Psi(0,\tau) ');
savefig('case2rate');
figure(5);
plot(x,lambda,'ro-',x,ones(1,N+1)*lambda_u,'k','linewidth',2,'markerfacecolor','r','markersize',6)
xlabel('\xi');
ylabel(' \lambda(\xi) ');
lamx =[' \lambda_u = ',num2str(lambda_u)];
legend(' Classical P3D',lamx)
if lambda(1)<lambda_u
    ax=axis;ax(3)=1;
    ax(4)=1.2*lambda(1);
    axis(ax);
end
pause(0.1)
end

function Omega = exact_PKN(x,t)
    gamma_0 = 1.0006328466775270;
    Omega_0 = ((12/5)*gamma_0^2)^(1/3);
    Omega   = Omega_0*t^(1/5)*(1-x).^(1/3).*(1-(1-x)/96 +23*(1-x).^2/64512);
end

function Omega = exact_PKNC(x,t,C)
    gamma_m0 = 2/pi/C; %#ok<NASGU>
    Omega_m0 = 2^(11/8)/pi^(1/2)/(3*C)^(1/4);
    Omega    = Omega_m0*t^(1/8)*(1-x).^(3/8).*(1+(1-x)/80);
end

[Walter Roberson]

Why did you get rid of the safety measures I put in for you in https://www.mathworks.com/matlabcentral/answers/596452-change-of-boundary-condition#comment_1027729 ?

Answer 1

[KSSV]

This line:

plot(t,Psi(0,t),'bo--','linewidth',2,'markerfacecolor','b')

t is a scalar value first i.e. t = 1.1000e-06; and Psi is an array of size 1X31....it looks like you are trying to access some value from Psi using: Psi(0,t), this is not correct; index cannot be zero, negative and fraction in MATLAB. You need to rethink on your code. Your code is buggy and got lot ov variables which are not being used.

Learn about debugging this will help you a a lot.

[Walter Roberson]

Other than changing some initialization constants such as number of iterations, and getting rid of some of the safety layers that I carefully added to ensure that the results remained reasonable, the only substantive difference I see is that you changed

Psi0 = (1-x).^(1/3);

to

 Psi0 = t0^(1/9)*(1-x).^(1/3);

In particular, you do not create any Psi function and you do not create anything that is obviously Psi(0,something) .

The program has an existing Psi0 variable; perhaps you wanted to plot that?

[Hai Nguyen]

Sorry Sir I happened to use the version i removed the safety layers to do some analysis while trying to understand the lines.

I want to plot time-dependent Psi profile at x=0. From what I understand Psi0 should be equal to Psi(0,t) but I'm unsure if t0=t from the whole codes. You can help me to explain this?

Does plot Psi0(0,t0) equal to Psi0(0,t)?

Is the following line for Psi computation?

Psi=[y(2,:) 0]

[Walter Roberson]

Psi(1) appears to be the value for x = 0, if I read the code correctly.

plot(t,Psi(1),'bo--','linewidth',2,'markerfacecolor','b')

[Hai Nguyen]

The resulting chart using the above line is below. It should be a curve to illustrate Psi profile in time at x = 0.

[Walter Roberson]

No it should not -- you have no "hold on" and t is scalar, so you are only plotting a single point and removing the one point the next time you go to plot in the same axes.

[Hai Nguyen]

Can you please explain why the points are not connecting to each other although i used 'o--'?

[KSSV]

It seems you are plotting point by point, like shown below:

figure
hold on
for i = 1:10
    plot(i,rand,'.')
end

To join them, don't plot point by point. Save the points into an array and then plot them at once. Like below:

a = zeros(10,1) ; 
for i = 1:10
a(i) = rand ;
end
plot(1:10,a)