This is interesting, but I am afraid it is too complicated to grasp in one sitting. One thing I can say is that sometimes I have been able to improve performance by many orders of magnitude by thinking about the problem differently. I mean, thinking of a different way to ask Mathematica the same question. But that is just generality. The only specific I can give right now is, as a first cut, to try using DiscretePlot:

DiscretePlot[f[x], {x, xmin, xmax, dx}]

This saves Mathematica the time to make decisions about where to evaluate, and it can give you a quick look at a complicated plot faster than you would get using Plot. Make sure you set Filling to False, otherwise the graph will look different than you expect.

I have a background in physics. If you send a reference to what you are trying to do written in human language, I may be able to digest it better. Best,

OL.

I sped it up another factor 10 just now, but there is still lots to go. I'm sure.

Attachments:

When you make the definitions, like you do know, it will evaluate it already and thus 'copying/replacing' the definition that occur in it. Since your definitions are very much nested, it will nest all these functions inside each other symbolically:

Compare the two ways of defining functions (immediate and delayed assignment):

ClearAll[a,b,n]
b[n_]:=n^2
a[m_]:=m b[m]
DownValues[a]
DownValues[b]

ClearAll[a,b,n]
b[n_]=n^2
a[m_]=m b[m]
DownValues[a]
DownValues[b]

Have you tried to evaluate T for some suitable numbers when you have all your definitions done using := (this way the kernel will not die of memory requirements because it does not substitute all the sub-definitions inside each other)

If you define everything with := rather than =, and just evaluate it numerically (machine numbers), it should be quite fast right? Doing this symbolically will probably be a nightmare anyhow, I quickly inspected the code, if you substitute everything inside each other it will be a HUGE expression. And evaluating this huge expression numerically probably takes longer than just evaluating it numerically directly. Have you tried?

I would not expect any sensical answer to appear with all these nested definitions.

I see only function definitions, what do you want to do with it in the end?

Without posting more details I'm afraid there is little that can be helped. We don't know SMz and Sz and L functions...

Doing MatrixPowers on symbolic expressions creates crazy big and complicated matrices very quickly, i think 'powering' numbers (not symbols) makes way more sense...

Please post the notebook!

Thanks for replying. I don't know how far down the thread you have read.

It's the basic square potential well problem. My tweaks are that I want a generic wave including spin, arbitrarily many identical barriers uniformly spaced ("a" apart and "d" wide) in graphene, with spin orbit interaction and a zeeman field under the potential.

My approach is to first try to replicate some simplified systems (only spin orbit, only one barrier, normal incidence etc) and then to generalise entirely so I can have any number of barriers at any incidence angle and with both SOI and zeeman.

I understand the maths of it (for the most part) but I am having huge trouble translating it to mathematica language cos I am far from fluent. I am making all of the rookie mistakes and it's costing me computing time, memory, and accuracy!

Delayed evaluation works for me, but only if you delay everything but the function that rely on Eigensystem commands...

<<Progress Update>>

Now I have successfully expanded this for two barriers (or a square potential well) and am attempting to expand fro arbitrarily many barriers. Getting errors I don't understand though. Also still can't delay evaluation without getting the wrong result.

Attachments:

I'm not quite sure, but I guess generally it has to do with cases there the numbers become '0' and you have to do some division. So is it possible for your elements in your matrix to be zero? Just check what you expect, you can put some values in there and try. I put 3 random values for qx qy and l in there...

Also you can check the symbolic result, is that what you want?

Have a look at the following functions:

AbsArg
AngleVector
ReIm
ArcTan[x, y]

The first gives the magnitude and the 'angle' (argument) of a complex number at once. The second can construct a vector by a radius-angle pair The third gives you real and Imaginary part of a number as a pair. And perhaps use Arctan with two arguments; so it takes care of the 4 quadrants (the signs of x and y et cetera).

The problem is in your eigensystem function:

First symbolic finding the general eigenvalues and eigenvectors may involve certain assumptions, not valid for all cases! Take for example the simple case of this polynomial:

In the general case 3x^2/x is 3x but not for x=0 ! So 'underneath' the symbolic calculation there are always some assumptions (unavoidable in some sense). I'm not sure why the general solution of the eigensystem 'works' for your input but not the non-general case!

Typical values I'm trying now

n=5. (this really should be exactly 5 but if this runs faster 5.0 will do)

a=0.23

d=0.12

kx= 1. (this is what I want to plot over so 1. is fine for now)

l=0.01

g=0.09