@testht06 What do you mean by "don't run"!? Do you mean that it doesn't stop running? Does it give an error message? Does simply stop running, but without producing any answer?

@testht06 What exactly happened when you ran it? You said that the program "don't run." For me, it doesn't stop running. I've let the Primfield command run for five hours of CPU time, and there's no result. Is this what happened to you?

I don't know if it's a bug in Primfield or if the computation really takes a long time.

For the other readers: It is computing the splitting field of an irreducible degree-three polynomial over GF(2^8). So the final field is only degree 24. Could it legitimately take more than five hours to compute it? The specific command is

alias(x= RootOf(Z^8+Z^7+Z^6+Z+1)):
alias(
     q= RootOf(
          Z^3 + (x^7+x^6+x^4+1)*Z^2 + (x^3+x+1)*Z + (x^7+x^6+x^4+x^3+x^2+x+1)
     )
):
Primfield({x,q}) mod 2;

The memory usage is low and stable: It's been at 66 Mb since the first few minutes.

@roman_pearce It's certainly more elegant than my solution. Vote up. Since the command inner is undocumented (see ?UndocumentedNames) and built-in, would you please explain what it does?

I think that your procedure would also benefit by doing a preliminary Expand.

Your code looks like Matlab.

Any subsindets of the form

subsindets(A, identical(B), C)

is equivalent to

subs([B]=~ C~([B]), A)

So, your inner subsindets could be simplified to

subs(alpha= freeze(alpha), x)

@acer Agreed: It would be better to create the final Matrix straightaway. It looks like the OP is trying to create a banded matrix, so the special input techniques for those would be helpful.

@billyp245 With powers of the variables a and b being arguments of erf and BesselK, these equations are far, far too complicated to solve symbolically. Your only hope is for a numeric solution, but for that you'll need to give numeric values to the constants.

@Art Kalb I think that they are equivalent for polynomial input. You should include the call to Expand within your procedure mod2simp. Both techniques will give incorrect results for unexpanded polynomials. For example, (1+alpha)^2 would become 1+alpha.

@Art Kalb Sorry that I didn't spot that. The whole purpose of the inner subsindets is to protect those alphas (by freezing them).

This works:

thaw(
     subsindets[2](
          subsindets[flat](z2, {identical(alpha), identical(alpha)^anything}, freeze),
          `^`, 1, op
     )
);

To be a bit more general, we can freeze any RootOf. Just replace identical(alpha) with specfunc(RootOf):

thaw(
     subsindets[2](
          subsindets[flat](z2, {specfunc(RootOf), specfunc(RootOf)^anything}, freeze),
          `^`, 1, op
     )
);

@billyp245 What do you mean by both derivatives? The expression D2(a,b) can only be one derivative. If you had two derivatives, say D2__a and D2__b, then you could solve for them being 0 with

solve({D2__a, D2__b}(a,b), {a,b});

If you have another question, please upload a worksheet.

For what it's worth, continuation doesn't help with this problem. I was able to successfully apply continuation (I multiplied the right side of the first boundary condition by the continuation parameter) only to get the much-more-difficult-to-deal-with error Newton iteration is not converging.

I need to know the modulus, i.e., the degree-8 irreducible polynomial, that you used to create GF(2^8).

Here's Kitonum's procedure with my improvements.

Download ContoursWithLabels.mw

@Kitonum Okay, I made it an Answer, and I'll make it a Comment to your post.

I did just what I said: changed x to x(t) and y to y(t). Then you just press ENTER and you should get the plot.