Scilab.tex

\documentclass[12pt]{report}
\usepackage{graphicx}
\begin{center}
 \huge{Computer Simulation in Physics}\\
 \huge{PAS 414}
\end{center}

\begin{figure}[h!]
  \centering
 \includegraphics[width =10cm,angle=360]{enis8.eps}
\end{figure}
\begin{table}[h!]
\centering
 Name : Lalit Kumar\\
 Roll No. : {\bf{CUHP17PGPAS11}}\\
 Submitted to : \bf{Dr. Jyoti Bhardwaj}
\end{table}
~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\\
\\

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{\underline{{\bf{Signature}}}}...................
\tableofcontents
\listoffigures
\begin{document}
 \chapter{Scilab Basics}
 \section{Basic Plots}
 Here I am ploting sine function in interval [0,1] with 10 sample points. Below is the Scilab
 commands that I used.\\
 \begin{verbatim}
  
-->a=0;b=1;N=10;
 
-->x = linspace(a,b,N);
 
-->y = sin(x);
 
-->plot(x,y,"k")
 
 \end{verbatim}

 \begin{figure}[h!][ht]
  \centering
\includegraphics[width=10cm,,angle=360]{enis2.eps} 
  \caption{Plot of Sin(x)}
  \end{figure}
  Now I am plotting this in interval [0,{$2\pi$}] with same points and different sample points 
  i.e., 30 points. Following commands I used\\
  \begin{verbatim}
   
-->a=0;b=2*%pi;N=10;
 
-->x = linspace(a,b,N);
 
-->y = sin(x);
 
-->plot(x,y,"k")
 
-->a=0;b=2*%pi;N=30;
 
-->x = linspace(a,b,N);
 
-->y = sin(x);
 
-->plot(x,y,"k")
 
  \end{verbatim}
Here is the plots.
\begin{figure}[h!]
 \centering
 \includegraphics[width = 10cm, ,angle=360]{enis1.eps}
 \caption{Plot of Sin(x) in interval [0,{2$\pi$}]} with 10 points.
\end{figure}

\begin{figure}[h!]
 \centering
 \includegraphics[width = 10cm, ,angle=360]{enis.eps}
 \caption{Plot of Sin(x) in interval [0,{2$\pi$}]} with 30 points.
\end{figure}
{\bf{Exercise 1.} Generate a variable s from [0,1]. Determine the values of the function 
{$y=sin{2\pi n s}$} for n = 0.1 and 10. Plot the functions.}\\
{\bf{Solution.}} Type following commands \\
\begin{verbatim}
 ->a=0;b=1;N=100;
 
-->s = linspace(a,b,N);
 
-->n=0.1;
 
-->y = sin(2*%pi*n*s);

-->plot(s,y,"k")
 
-->n=10;
 
-->a=0;b=1;N=100;
 
-->s = linspace(a,b,N);
 
-->y = sin(2*%pi*n*s);
 
-->plot(s,y,"k")
\end{verbatim}
Output from these commands is
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{enis3.eps}
 \caption{Plot of {$y=sin(2\pi n s)$} at n= 0.1}
\end{figure}

\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{enis4.eps}
 \caption{Plot of {$y=sin(2\pi n s)$} at n= 10}
\end{figure}
{\bf{Exercise 2.} Plot a sine wave with frequency 20Hz and initial phase {$\frac{\pi}{2}$} in the 
interval [0,1]}\\
{\bf{Solution}} Type following commands\\
\begin{verbatim}
 
-->a=0;b=1;N=100;
 
-->x = linspace(a,b,N);
 
-->f = 20;
 
-->phi = %pi/2;
 
-->y = sin(2*%pi*f*x+phi);

-->plot(x,y,"k")
\end{verbatim}
Graph of output is
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{enis5.eps}
 \caption{Plot of sine wave with f= 20 Hz,{$\phi = \frac{\pi}{2}$}}
\end{figure}
\section{To Study Interference of Waves}
Type following commands
\begin{verbatim}
 ->x = 0:0.01:1;
 
-->f = 8;
 
-->y1 = sin(2*%pi*f*x);
 
-->phi = 0;

 \\**************INITIAL PHASE = 0 ****************************
 
-->y2 = sin(2*%pi*f*x+phi);
 
-->subplot(2,2,1)
 
-->plot(x,y1+y2,"k")

\\ ****************FIRST PLOT *********************************

-->phi = %pi;

 \\*************************PHASE DIFF. = PI *******************
 
-->y3= sin(2*%pi*f*x+phi);
 
-->subplot(2,2,2)
 
-->plot(x,y1+y3,"k")

\\ *****************SECOND PLOT**********************************

-->phi = %pi/2;
 
 \\***********************PHASE DIFF. = PI/2 ********************
 
-->y4= sin(2*%pi*f*x+phi);
 
-->subplot(2,2,3)
 
-->plot(x,y1+y4,"k")
 
 \\**********************THIRD PLOT*******************************
 
-->phi = %pi/4;
 
 \\**********************PHASE DIFF. = PI/4 **********************
 
-->y5= sin(2*%pi*f*x+phi);
 
-->subplot(2,2,4)
 
-->plot(x,y1+y5,"k")
 
 \\*****************************FOURTH PLOT**************************
\end{verbatim}
Output is
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{enis6.eps}
 \caption{Interference}
\end{figure}
\section{Study of Beats}
Type following commands in Scilab
\begin{verbatim}
-->x = 0:0.001:1;
  
-->f1 = 19;
 
-->y1 = sin(2*%pi*f1*x);
 
-->f2 = 23;
 
-->y2 = sin(2*%pi*f2*x);
 
-->subplot(3,1,1)
 
-->plot(x,y1,"k")
 
-->subplot(3,1,2)
 
-->plot(x,y2,"k")
 
-->subplot(3,1,3)
 
-->plot(x,y1+y2,"k")
 
\end{verbatim}
\begin{verbatim}
 exec('/home/lalit/beats.sci', -1)
 
-->beats(0,1,1000);
 
-->A = locate(6)
 A  =
 
 
         column 1 to 5
 
    0.1786885    0.1896175    0.2355191    0.2879781    0.3732240  
    14.160104    14.160104    14.160104    14.160104    13.923219  
 
         column 6
 
    0.4191257  
    13.923219  
 
-->tavg = diff(A(1,:))
 tavg  =
 
    0.0109290    0.0459016    0.0524590    0.0852459    0.0459016  
 
-->favg = 1/mean(tavg)
 favg  =
 
    20.795455  
 
-->B = locate(4)
 B  =
 
  - 0.0857923    0.2158470    0.2377049    0.4825137  
  - 1.5216926    0.5155197    0.4207657    0.3733886  
 
-->tbeat = diff(B(1,:))
 tbeat  =
 
    0.3016393    0.0218579    0.2448087  
 
-->fbeat = 1/mean(tbeat)
 fbeat  =
 
    5.2788462  
 
\end{verbatim}

\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{enis7.eps}
 \caption{Beats}
\end{figure}
{\bf{Exercise : Modulation: Generate wave of frequencies where one of them is of very high
frequency(carrier) and the other with low frequency(signal). For example, choose {$f1 = 100 Hz$}
and {$f2 = 3 Hz$}. Determine the product of the two wavs y and y2 as y1.*y2(.* is the operation 
for obtaining element by element multiplication, while * would give the error as it is refered for matrix multiplication
in Scilab) }}\\
{\bf{Solution}} Type following commands in Scilab
\begin{verbatim}
 ->x = 0:.0001:1;

-->f1  = 100;
 
-->y1 = sin(2*%pi*f1*x);
 
-->f2 = 3;
 
-->y2 = sin(2*%pi*f2*x);

-->z = y1.*y2;
 
-->plot(x,z,'k')
 
 \end{verbatim}
 \begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{enis8.eps}
 \caption{Modulation}
\end{figure}
\section{Fourier Superposition Theorem}
Any periodic waveform {$f(x)$ } with periodicity L can be represented as linear combination 
of simple harmonic wave and their harmonics of appropriate amplitude and is mathematically
represented as:\\
\begin{equation}
 f(x) = \sum_{n=0}^{\infty} a_{n}cos{(\frac{2n\pi x}{L})} + b_{n}sin{(\frac{2n\pi x}{L})}
\end{equation}
where the Fourier coefficients {$a_{n}$} and {$b_{n}$} are given by\\
\begin{equationarray}
 a_{n} = \frac{1}{L}\int_{\frac{-L}{2}}^{\frac{L}{2}} f(x)cos{(\frac{2n\pi x}{L})}dx\\
 b_{n} = \frac{1}{L}\int_{0}^{2\pi} f(x)sin{(\frac{2n\pi x}{L})}dx
 
\end{equationarray}
\subsection{Fourier Series for a Square Wave}
A square wave with periodicity L is given by\\
\begin{equation}
 Sqr(x) = \sum_{n=1,3,5}^{\infty} b_{n}sin{(\frac{2n\pi x}{L})}
 \end{equation}
 with {$b_{n} = \frac{1}{n}$}. Here {$a_{n}$}s are all small.\\
 \section{Fourier Superposition}
 Type following commands in Scilab
 \begin{verbatim}
  -->x = 0:0.001:1;
 
-->f = 1;
 
-->y1 = sin(2*%pi*f*x);
 
-->f = 3 ;
 
-->y2 = sin(2*%pi*f*x);
 
-->f = 5 ;
 
-->y3 = sin(2*%pi*f*x);
 
-->y = y1 + y2/3 + y3/5 ;
 
-->subplot(2,2,1)
 
-->plot(x,y1)
 
-->subplot(2,2,2)
 
-->plot(x,y2)
 
-->subplot(2,2,3)
 
-->plot(x,y3)
 
-->subplot(2,2,4)
 
-->plot(x,y,'r')
 \end{verbatim}
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{Fsuper.eps}
 \caption{Plots of Square wave.}
\end{figure}
{\bf{Exercise:Add more number of components(say 5,7, and 10) to obtain better approximations of the square wave which will
be reflected as decrease in the amplitude of oscillation in the flat region as well as in slope
at the falling edges. Tabulate the results.}}\\
{\bf{Solution}} Type the following commands in Scilab
\begin{verbatim}
 -->x = 0:0.001:1;
 
-->f = 1;
 
-->y1 = sin(2*%pi*f*x);
 
-->f = 3;
 
-->y2 = sin(2*%pi*f*x);
 
-->f = 5;
 
-->y3 = sin(2*%pi*f*x);
 
-->f = 7;
 
-->y4 = sin(2*%pi*f*x);
 
-->f = 9;
 
-->y5 = sin(2*%pi*f*x);
 
-->f = 11;
 
-->y6 = sin(2*%pi*f*x);
 
-->f = 13;
 
-->y7 = sin(2*%pi*f*x);
 
-->f = 15;
 
-->y8 = sin(2*%pi*f*x);
 
-->f = 17;
 
-->y9 = sin(2*%pi*f*x);
 
-->f = 19;
 
-->y10 = sin(2*%pi*f*x);
 
-->y = y1 + y2/3 + y3/5 + y4/7 + y5/9 + y6/11 + y7/13 + y8/15 + y9/17 +y10/19;
 
-->subplot(2,4,1)
 
-->plot(x,y1)
 
-->subplot(2,4,1)
 
-->plot(x,y4)
 
-->subplot(2,4,2)
 
-->plot(x,y5)
 
-->subplot(2,4,3)
 
-->plot(x,y6)
 
-->subplot(2,4,4)
 
-->plot(x,y7)
 
-->subplot(2,4,5)
 
-->plot(x,y8)
  
-->subplot(2,4,6)
 
-->plot(x,y9)
 
-->subplot(2,4,7)
 
-->plot(x,y10)
 
-->subplot(2,4,8)
 
-->plot(x,y)
 
\end{verbatim}

\begin{figure}[h!]
 \centering
 \includegraphics[width =15cm,hight=15cm,angle=360]{Fsupermore.eps}
 \caption{Adding first 10 components}
\end{figure}
{\bf{Exercise : Repeat the above procedure to produce a square wave of frequency f = 2.5Hz}}\\
{\bf{Solution}} Type the following commands in Scialab
\begin{verbatim}
 ->x = 0:0.001:1;
 
-->f = 2.5;
 
-->y1 = sin(2*%pi*f*x);
 
-->f = 7.5;
 
-->y2 = sin(2*%pi*f*x);
 
-->f = 12.5;
 
-->y3 = sin(2*%pi*f*x);
 
-->y = y1 + y2/3 + y3/5;
 
-->subplot(2,2,1)
 
-->plot(x,y1)
 
-->subplot(2,2,2)
 
-->plot(x,y2)
 
-->subplot(2,2,3)
 
-->plot(x,y3)
 
-->subplot(2,2,4)
 
-->plot(x,y)
 
\end{verbatim}
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{Sq25hz.eps}
 \caption{Plots of Square wave of frequency 2.5Hz}
\end{figure}
\section{Program to add n terms in the Fourier Series}
Open SciNotes using commands
\begin{verbatim}
 ->edit SquareWave.sci
\end{verbatim}
\begin{figure}[h!]
 \centering
 \includegraphics[width =12cm,hight=17cm,angle=360]{Sq25hz10tr1.eps}
 \caption{Command in SciNotes.}
\end{figure}
Output is
\\
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{Sq25hz10tr.eps}
 \caption{A square wave of 2.5 Hz obtained by adding first 10 terms in its Fourier Series expansion}
\end{figure}
\chapter{Wave Packet}
\section{Construction of a Wave-packet using Superposition}
Type this command to open a SciNotes\\
\begin{verbatim}
 
-->edit wavepacket.sci
 
\end{verbatim}
And type commands 

\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{wavepacketsci.eps}
 \caption{Commands in SciNotes}
\end{figure}
Output is as
\begin{verbatim}
 ->exec('/home/lalit/wavepacket.sci', -1)
 
-->y1 = wavepacket(1,10,1)
\end{verbatim}

\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{wavepacket.eps}
 \caption{Superposition of waves: spacing b/w the waves is dfs = 1 Hz, number of waves added is n = 10 and the 
 region is x = [-1,1]}
\end{figure}
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,angle=360]{wavepacketpar.eps}
 \caption{Wave-packet construction for different parameters.}
\end{figure}
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,angle=360]{wavepackettyp.eps}
 \caption{A typical wavepacket}
\end{figure}
\subsection{Gaussian amplitude function}
The Gaussian function is given by
\begin{equation}
 a(k) = \frac{1}{\sqrt{(2\pi (\Delta k)^2)}}exp{(-\frac{(k-k_{0})^2}{2\Delta k^2})}
\end{equation}
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,angle=360]{Gaussian.eps}
 \caption{Wavepacket for a Gaussian amplitude function}
\end{figure}
\subsection{2-D Gaussian Wave-packet}
\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,angle=360]{Gaussian2D.eps}
 \caption{Wavepacket for a 2D Gaussian amplitude function}
\end{figure}

\section{Sawtooth and Triangular waves }
\begin{figure}[h!]
 \centering
 \includegraphics[width = 10cm,angle = 360]{Sawtooth_wave.eps}
 \caption{Plot of Sawtooth wave}
\end{figure}

\begin{figure}[h!]
 \centering
 \includegraphics[width = 10cm,angle = 360]{Trianguler_wave.eps}
 \caption{Plot of Triangular wave}
\end{figure}

\begin{figure}[h!]
 \centering
 \includegraphics[width = 10cm,angle = 360]{Sawtooth_waveScript.eps}
 \caption{SciNotes commandsfor Sawtooth Wave}
\end{figure}

\begin{figure}[h!]
 \centering
 \includegraphics[width = 10cm,angle = 360]{Triangular_waveScript.eps}
 \caption{SciNotes Commands for Triangular Wave}
\end{figure}
\chapter{Matrixmethod to solve Schrodinger equation for different kind of Potentials}
\section{Introduction} A well known technique to find numerical solution of Schrodinger equation i.e., energy eigenvalues//
 is done by the implication of Matrixmethod. To reach at these solutions Scilab is choosen as programming language. \\
 
 \section{Scilab code for a Squarewell embedded in infinite squarewell}
 \begin{verbatim}
  function [E,h]=squarewell(V0,b,a,N,q)
    a = a/0.52917725;
    b = b/0.52917725;
    V0= V0/27.211396;
    x = 0:0.001:a;
    n=length(x);
    V=V0*ones(1,n);
    b2=n/2;
    b1=(n/2)-((n*b)/(2*a));
    b3=(n/2)+((n*b)/(2*a));
    V(b1:b3)=0;
    mprintf('b1=%f,b2=%f,b3=%f',b1,b2,b3);
    h=zeros(N,N); //initialising the h matrix 
    s=eye(N,N);   //equivalent identity matrix
    for i=1:N
        for j = i:N
            bi = sqrt(2/a)*sin(i*%pi*x/a);
            bj = sqrt(2/a)*sin(j*%pi*x/a);
            f = bi.*V.*bj;
            h(j,i) = intsplin(x,f);
            if(j==i) then
                h(j,i)=h(j,i)+j^2*%pi^2*(0.5/a^2);
            end 
            h(i,j)=h(j,i);
        end
    end
    [al,bl,R]=spec(h,s); //Solving eigen value problem
    e1=al./bl;           //Eigen values
    c=R;                 //Eigen vectors
    [E,k]=gsort(e1);    //Sorting the eigen values
    E=E*27.211396;
    d=c;
    d(:,1:N)=c(:,k);    //Corresonding eigen vectors
    j=N;//variable to keep track of wavefn
    i=1;//loop variable for generating jth wavefn
    while j>=1
        psi(j,1:n) = 0; //initialising the wavefn
        m=1;
        while m<=N
            psi(j,:)=psi(j,:)+d(m,i)*sqrt(2/a)*sin(m*%pi*x/a);//wavefn using eigen expansion
            m=m+1;
        end
        i=i+1;   
        j=j-1;
     end 
     p=1;//loop variable for choosing the pth wavefn to be plotted
     x=x*0.52917725;
     plot(x,V*27.211396)
     while p<=q*2
         plot(x,psi(p,:)+E(N-p+1)*ones(1,n));
         plot(x,E(N-p+1)*ones(1,n),"r")
       // xclick();
         //clf()
         p=p+1;
     end
endfunction
 \end{verbatim}
 \begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{foo1.eps}
 \caption{Wavefunctions}
\end{figure}

\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,angle=360]{foo1p.eps}
 \caption{Probability}
\end{figure}

 \section{Scilab code for a Double Squarewell embedded in infinite squarewell}
 \begin{verbatim}
  function [E,h]=matrixdoublesquarewell(V0,b,a,c,N,q)
    //Vo is depth of the well
    //b is width of the square well
    //a is width of the infinite 1D potential 
    //N is the number of terms in the basis fn expansion
    //q is the number of wavefns to be plotted
    a = a/0.52917725;
    b = b/0.52917725;
    c=c/0.52917725;
    V0= V0/27.211396;
    x = 0:0.001:a;
    n=length(x);
    V=V0.*ones(1,n);
    b3=n/2;
    b1=(n/2)-((n*c)/(2*a))-((b*n)/a);
    b2=(n/2)-((n*c)/(2*a));
    b4=(n/2)+((n*c)/(2*a));
    b5=(n/2)+((n*c)/(2*a))+((b*n)/a);
    V(b1:b2)=0;
    V(b4:b5)=0;
    mprintf('b1=%f,b2=%f,b3=%f,b4=%f,b5=%f',b1,b2,b3,b4,b5);
 
    h=zeros(N,N); //initialising the h matrix 
    s=eye(N,N);   //equivalent identity matrix
        
    //Determing the elements of the h matrix
    for i=1:N
        
        for j = i:N
            bi = sqrt(2/a)*sin(i*%pi*x/a);
            bj = sqrt(2/a)*sin(j*%pi*x/a);
            f = bi.*V.*bj;
            h(j,i) = intsplin(x,f);
            if(j==i) then
                h(j,i)=h(j,i)+j^2*%pi^2*(0.5/a^2);
            end 
            h(i,j)=h(j,i);
        end
    end
    [al,bl,R]=spec(h,s); //Solving eigen value problem
    e1=al./bl;           //Eigen values
    c=R;                 //Eigen vectors
    [E,k]=gsort(e1);    //Sorting the eigen values
    E=E*27.211396
    d=c;
    d(:,1:N)=c(:,k);    //Corresonding eigen vectors
    
    j=N;//variable to keep track of wavefn
    i=1;//loop variable for generating jth wavefn
    while j>=1
        psi(j,1:n) = 0; //initialising the wavefn
        m=1;
        while m<=N
            psi(j,:)=psi(j,:)+d(m,i)*sqrt(2/a)*sin(m*%pi*x/a);//wavefn using eigen expansion
            m=m+1;
        end
        i=i+1;   
        j=j-1;
     end 
     p=1;//loop variable for choosing the pth wavefn to be plotted
     x=x*0.52917725;
     plot(x,V*27.211396)
     while p<=q
         plot(x,psi(p,:)+E(N-p+1)*ones(1,n));
         plot(x,E(N-p+1)*ones(1,n),"r")
        // xclick();
        // clf()
        
         p=p+1;
     end
endfunction
 \end{verbatim}

\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,angle=360]{foo1pp.eps}
 \caption{Wavefunctions}
\end{figure}

\begin{figure}[h!]
 \centering
 \includegraphics[width =10cm,,angle=360]{foo1pr.eps}
 \caption{Probability}
\end{figure}

 \section{Scilab code for a n-Squarewells embedded in infinite squarewell}
 \begin{verbatim}
  function [E,h]=unsquarewell(V0,d,b,a,N,q)
    a = a/0.52917725;
    b = b/0.52917725;
    V0= V0/27.211396;
    x = 0:0.001:a;
    n=length(x);
    V=V0*ones(1,n);
    period = a/(d+1);
    for i = 1:d
        t = (i*period-b/2)*n/a;
        s = (i*period+b/2)*n/a;
        V(t:s) = 0;
    end
    plot(x*0.52917725,V*27.211396)
    h=zeros(N,N); //initialising the h matrix 
    s=eye(N,N);   //equivalent identity matrix
    for i=1:N
        for j = i:N
            bi = sqrt(2/a)*sin(i*%pi*x/a);
            bj = sqrt(2/a)*sin(j*%pi*x/a);
            f = bi.*V.*bj;
            h(j,i) = intsplin(x,f);
            if(j==i) then
                h(j,i)=h(j,i)+j^2*%pi^2*(0.5/a^2);
            end 
            h(i,j)=h(j,i);
        end
    end
    [al,bl,R]=spec(h,s); //Solving eigen value problem
    e1=al./bl;           //Eigen values
    [E,k]=gsort(e1);    //Sorting the eigen values
    E=E*27.211396;
    d(:,1:N)=R(:,k);    //Corresonding eigen vectors
    
    j=N;//variable to keep track of wavefn
    i=1;//loop variable for generating jth wavefn
    while j>=1
        psi(j,1:n) = 0; //initialising the wavefn
        m=1;
        while m<=N
            psi(j,:)=psi(j,:)+d(m,i)*sqrt(2/a)*sin(m*%pi*x/a);//wavefn using eigen expansion
            m=m+1;
        end
        i=i+1;   
        j=j-1;
     end 
     p=1;//loop variable for choosing the pth wavefn to be plotted
     x=x*0.52917725;
   
     while p<=q
         plot(x,psi(p,:)+E(N-p+1)*ones(1,n));
         plot(x,E(N-p+1)*ones(1,n),"r")
       // xclick();
         //clf()
        
         p=p+1;
     end
endfunction
 \end{verbatim}

  \section{Scilab code for a harmonic Oscillator embedded in infinite squarewell}
  \begin{verbatim}
   function [E,h]=harmonic(k,a,N,q)
	//Vo is depth of the well
	//b is width of the square well
	//a is width of the infinite 1D potential 
	//N is the number of terms in the basis fn expansion
	//q is the number of wavefns to be plotted
	a = a/0.52917725;
   	// b = b/0.52917725;
    	//c = c/0.52917725;
   	// V0= V0/27.211396;
    	k = k/0.52917725;
    	x = 0:0.001:a;
    	n=length(x);
    	V=zeros(1,n);
   	//mprintf('b1=%f,b2=%f',b1,b2);
    	V=(1/2)*k*(x-a/2).*(x-a/2);  
    	h=zeros(N,N); //initialising the h matrix 
    	s=eye(N,N);   //equivalent identity matrix
        
    	//Determing the elements of the h matrix
    	for i=1:N
        
        	for j = i:N
            	bi = sqrt(2/a)*sin(i*%pi*x/a);
            	bj = sqrt(2/a)*sin(j*%pi*x/a);
            	f = bi.*V.*bj;
            	h(j,i) = intsplin(x,f);
            	if(j==i) then
                	h(j,i)=h(j,i)+j^2*%pi^2*(0.5/a^2);
            	end 
            	h(i,j)=h(j,i);
        	end
    	end
    	[al,bl,R]=spec(h,s); //Solving eigen value problem
    	e1=al./bl;           //Eigen values
    	c=R;                 //Eigen vectors
    	[E,k]=gsort(e1);    //Sorting the eigen values
    	E=E*27.211396;
    	d=c;
    	d(:,1:N)=c(:,k);    //Corresonding eigen vectors
    
   	 j=N;//variable to keep track of wavefn
    	i=1;//loop variable for generating jth wavefn
    	while j>=1
        	psi(j,1:n) = 0; //initialising the wavefn
        	m=1;
        	while m<=N
            	psi(j,:)=psi(j,:)+d(m,i)*sqrt(2/a)*sin(m*%pi*x/a);//wavefn using eigen expansion
           	 m=m+1;
        	end
       		i=i+1;   
        	j=j-1;
     	end
     	subplot(2,1,1)
      	p=1;//loop variable for choosing the pth wavefn to be plotted
     	x=x*0.52917725;
     	plot(x,V*27.211396)
    
     	while p<=q
         	subplot(2,1,2)
         	plot(x,psi(p,:)+E(N-p+1)*ones(1,n));
                plot(x,E(N-p+1)*ones(1,n),"r")
         	xclick();
         	clf();
         	p=p+1;
     	end
endfunction

  \end{verbatim}
  \section{Scilab code for Anharmonic Oscillator embedded in infinite squarewell}
  \begin{verbatim}
   function [E,h]=harcinom(k,alpha,a,N,q)

    a = a/0.52917725;
   // b = b/0.52917725;
    //c = c/0.52917725;
   // V0= V0/27.211396;
   alpha = alpha/0.52917725;
   k = k/0.52917725;
    x = 0:0.001:a;
    n=length(x);
    V=zeros(1,n);
  // mprintf('b1=%f,b2=%f',b1,b2);
    V=(1/2).*k.*(x-a/2).*(x-a/2)+alpha.*(x-a/2).^4;  
    h=zeros(N,N); //initialising the h matrix 
    s=eye(N,N);   //equivalent identity matrix
        
    //Determing the elements of the h matrix
    for i=1:N
        
        for j = i:N
            bi = sqrt(2/a)*sin(i*%pi*x/a);
            bj = sqrt(2/a)*sin(j*%pi*x/a);
            f = bi.*V.*bj;
            h(j,i) = intsplin(x,f);
            if(j==i) then
                h(j,i)=h(j,i)+j^2*%pi^2*(0.5/a^2);
            end 
            h(i,j)=h(j,i);
        end
    end
    [al,bl,R]=spec(h,s); //Solving eigen value problem
    e1=al./bl;           //Eigen values
    c=R;                 //Eigen vectors
    [E,k]=gsort(e1);    //Sorting the eigen values
    E=E*27.211396;
    d=c;
    d(:,1:N)=c(:,k);    //Corresonding eigen vectors
    
    j=N;//variable to keep track of wavefn
    i=1;//loop variable for generating jth wavefn
    while j>=1
        psi(j,1:n) = 0; //initialising the wavefn
        m=1;
        while m<=N
            psi(j,:)=psi(j,:)+d(m,i)*sqrt(2/a)*sin(m*%pi*x/a);//wavefn using eigen expansion
            m=m+1;
        end
        i=i+1;   
        j=j-1;
     end
     subplot(2,1,1)
      p=1;//loop variable for choosing the pth wavefn to be plotted
     x=x*0.52917725;
     plot(x,V*27.211396)
    
     while p<=q
         subplot(2,1,2)
         plot(x,psi(p,:)+E(N-p+1)*ones(1,n));
         plot(x,E(N-p+1)*ones(1,n),"r")
         xclick();
         clf();
         p=p+1;
     end
endfunction
  \end{verbatim}

\end{document}