# 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 = ...