# Octave is a free version of Matlab. 

# Functions like Psi0 are created with 
#    function y=...
# 	y = ...;
#    endfunction
# If you leave off the semicolon, it'll print out the value of the
# function each time it's called.

function amp = Psi0 (x)
  amp = (1/sqrt(3.)) * cos(... / (... *3.));
endfunction

Psi0(0.)

# Numerical integration from 'min' to 'max' is done using the function 'quad'.
# Its arguments are # the function name in quotes, followed by min and max. 
# Let's test whether we got the normalization right on our ground state,
# by integrating Psi0(x)^2 from -a to a, with a=3. 

# First, define the integrand:

function sq = Psi0Sq (x)
  sq = ... ^2;
endfunction 

Psi0Sq(0.)

# Then call quad. It returns not only the integrak, but some diagnostics
# and an estimated error:

[integral, ier, nfun, err] = quad ("Psi0Sq", -3, 3)

# We now define the trial wavefunction, without the right normalization
# Note: Many functions in Matlab like times have a version * for 
# multiplying numbers, and a version .* for termwise multiplication of 
# vectors. We'll use the .^ operation to square our position:

function amp = PsiUnnormalized (x)
   amp = 3.^2 - ...;
endfunction

# This is useful for plots and stuff, where you can define an array x 
# using 'linspace' and then plot Psi(x) with one function call, passing
# in the array. 

xCoarse = linspace(-3.,3.,7)
PsiUnnormalized(xCoarse)

function sq = PsiUnnormalizedSq (x)
  ...

[InverseNormSq, ier, nfun, err] = quad(...)
norm = ...

# Passing norm into the function Psi is apparently tricky in Matlab.
# We'll type in the value of norm explicitly...

function amp = Psi (x)
   amp = 0.0... * PsiUnnormalized(x);
endfunction

# Testing that the normalization is right...

function sq = PsiSq (x)
    ...
endfunction

[integral, ...

# Plotting the function at points x defined by linspace. Plot takes
# two arrays x and y, but since we've used termwise exponentiation in Psi(x)
# we define an array x and then pass it x, Psi(x).

# To plot two functions f(x) and g(x) on the same graph, use
# plot(x, f(x), x, g(x))

x = linspace(-3.,3.,100);
plot(x, Psi0(x), x, Psi(x))

# Define the integrand for the overlap of Psi0 and Psi:

function dot = OverlapIntegrand (x)
   ...
endfunction

[overlap, ier, nfun, err] = quad("OverlapIntegrand", -3, 3)

# Finally, find the probability that Psi is in the ground state

ProbGroundState = ...