% mesh
    DL = 0.08 ;
    sd_ratio = 3;
    SL = DL*sd_ratio ;
    Np = 8 ;
    N = floor(Np/DL);
    dx = L/N;
    a = dx ;
    x = zeros(1,N+1) ;
    q = ones(1, N) ;
    G = ones(1, N) ;
    lam = 1e-5 ; %
    beta = ro/EI ;
    dt= sqrt(lam*beta)*dx^2 ;

    % Boundaries
    yp1(N+1) = 0 ;
    y = 0 ;
    ym1 = 0 ;

    % Implementation
    i = 1 : N+1 ;
    x = dx * (i-1) ;

    tf = 0.1 * 1/f ;
    Nt = floor(tf/dt) ;
    for j=1:Nt
        t=dt*(j-1);
        for i=1:N
            if j==1, yp1(i)=0; %eq (7)
            elseif j==2
                if i==1
                    yp1(1) = EI/2/ro*dt^2/dx^4* ( 2*dx^3*p*cos(w*t)/EI ); %eq (8)
                elseif i==2,
                    G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
                    q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
                    yp1(i)= y(i) - q(i)/2/ro*dt^2 -EI/2/ro *dt^2/dx^4*(y(i+2) -4*y(i+1)+7*y(i)-4*y(i-1));
                    %eq(23)>> 6
                elseif i==N, yp1(i)=0;
                    %eq(25)>>10
                else
                    G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
                    if ( ( 0  <= rem( abs(x(i)),SL*x_f ) ) && ( rem(abs(x(i)),SL*x_f ) <= DL/2*x_f  ) ) || ( ( ( SL-DL/2)*x_f  <= rem( abs(x(i)),SL*x_f ) ) && ( rem( abs(x(i)),SL*x_f ) <= SL*x_f ) )
                        q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
                    else
                        q(i) = ks0*y(i)/( u*(1 + ks0*y(i)/qu )) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% Within the soil
                    end

                    yp1(i) = y(i) - q(i)/2/ro*dt^2 - EI/2/ro*dt^2/dx^4*(y(i+2)-4*y(i+1)+6*y(i)-4*y(i-1)+y(i-2)) ;
                    %eq(20)>>4
                end
            else
                if i==1
                    G(1) = g ;
                    q(1) = kc0*y(1)/( 1 + kc0*y(1)/qcu ) - G(1)*H*( 2*y(2) - 2*y(1) )/a^2 ;% within the sc

                    yp1(1) = 2*y(1) - ym1(1) - q(1) /ro*dt^2 - EI/ro * dt^2/dx^4*(2*y(3)-8*y(2) + 6*y(1)-2*dx^3*p/EI*cos(w*t) );
                    %eq (7)
                elseif i==2,
                    G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
                    q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc

                    yp1(i) = 2*y(i) - ym1(i) - q(i)/ro *dt^2 - EI/ro*dt^2/dx^4*(y(i+2) -4*y(i+1)+7*y(i)-4*y(i-1) );
                    %eq(22)>>5
                elseif i==N,
                    G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
                    if ( ( 0  <= rem( abs(x(i)),SL*x_f ) ) && ( rem(abs(x(i)),SL*x_f ) <= DL/2*x_f  ) ) || ( ( ( SL-DL/2)*x_f  <= rem( abs(x(i)),SL*x_f ) ) && ( rem( abs(x(i)),SL*x_f ) <= SL*x_f ) )
                        q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
                    else
                        q(i) = ks0*y(i)/( u*(1 + ks0*y(i)/qu )) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% Within the soil
                    end
                    yp1(i)= 2*y(i) - ym1(i) - q(i)/ro*dt^2 - EI/ro*dt^2/dx^4*(y(i-2) - 4*y(i-1)+5*y(i) );
                    %eq(24)>>9
                else
                    G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
                    if ( ( 0  <= rem( abs(x(i)),SL*x_f ) ) && ( rem(abs(x(i)),SL*x_f ) <= DL/2*x_f  ) ) || ( ( ( SL-DL/2)*x_f  <= rem( abs(x(i)),SL*x_f ) ) && ( rem( abs(x(i)),SL*x_f ) <= SL*x_f ) )
                        q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
                    else
                        q(i) = ks0*y(i)/( u*(1 + ks0*y(i)/qu )) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% Within the soil
                    end
                    yp1(i)=2*y(i)-ym1(i) - q(i)/ro*dt^2 - EI/ro*dt^2/dx^4*( y(i+2) -4*y(i+1)+6*y(i)-4*y(i-1)+y(i-2));
                    %eq(17)>>1
                end
            end
        end
        ym1 = y;
        y = yp1;

    end