The Euler method of solving differential equations is less accurate in predicting the solution than those used in the previous code. This required us to use a smaller time step of 100 seconds for time advancing.
First, I tried only the meridional flow without any diffusion or differential rotation disturbing the magnetic field. And to my surprise, I found, that not only there was no imbalance between the hemispheres, but even the flux within each hemisphere was conserved individually!
It is logically obvious, that for a symmetric initial condition, the flux distribution at any following time should be symmetric across the equator. And so it was!
Then I gave a try to the code with only differential rotation. Still no imbalance, and accurate flux conservation!
The next step was diffusion. The diffusion only code also conserved the flux individually in hemispheres, with almost no imbalance. Sometimes, there would be an imbalance, about 10 million times less than the flux in a sunspot. But then it would become zero again.
And then, I gave a run to the code, with all the 3 effects combined, and following is a screenshot showing:
1st column: Southern Hemisphere Flux
2nd column: Northern Hemisphere Flux
3rd column: Imbalance.
I also took the liberty of making a movie of the simulation, to see the flux reach the poles in 3D! The 3D plotting was one of the best things I did while in Montana.. :)
This movie was made from a simulation that runs for about 2 years in the code, and ran for about 24 hours on my laptop!