Nuclear Reactor Simple 1D Method of Characteristics Solver

[ChasingNeutrons]

Homework Statement

Hi all, attempting to make a Method of Characteristics solver in Matlab. I'm particularly hoping that there are some computational nuclear engineering guys about who might have a bit of experience programming a simple version of one of these.

I'm trying to create the solver to find the scalar neutron flux and criticality in a one-group, homogeneous, multiplying with slab subject to vacuum boundary conditions but, infuriatingly, while I can attain a neutron flux profile which is correct (verifying with other codes) my criticality is off by a factor of 10! I was hoping that someone may be able to provide a bit of insight into what I'm doing incorrectly - it may be a long shot but it's a problem I've been suffering with for far too long. If anyone would care to have a look at the code as it stands then do let me know!

Homework Equations

For the method of characteristics, the fundamental equations that I'm iterating around in my solver are quite standard. I won't include them here as I would rather just share my code if you are familiar with the general problem.

The Attempt at a Solution

For my problem, I'm examining a geometry of length 100cm with a Σ_tr of 0.2758/cm, νΣ_f of 0.0366/cm and Σ_a of 0.0253/cm. As I say, I can get the correct normalised neutron flux profile across the slab but sadly k is approximately a factor of 10 smaller than the k_eff value of about 1.3876 that I am expecting. It is interesting to note that for perfectly reflective boundary conditions I attain a k_∞ which is equal to νΣ_f/Σ_tr. As I say, feel free to ask for the code if you think you might be of any help - if you can identify the problem I'll owe you my dissertation!