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.

Don't repeat your Questions. It's very annoying to the moderators.

You got an Answer in about 35 minutes. I think that that's a pretty good response rate.

@WernerP There's no need for implicitplot3d or high values of numpoints. Your code can be simplified to

A:= plot3d(5, theta= -Pi..Pi, phi= arccos(4/5)..arccos(3/5), coords= spherical):
B:= plot3d(5, theta= -Pi..Pi, phi= 0..Pi/2, coords= spherical, transparency= .6):
C1:= plot3d([3, theta, z], theta= -Pi..Pi, z= 0..4, coords= cylindrical):
C2:= plot3d([4, theta, z], theta= -Pi..Pi, z= 0..3, coords= cylindrical):
plots:-display(
     {A, B, C1, C2}, scaling= constrained, labels= [X,Y,Z],
     axes= normal, style= patchnogrid, color= yellow
);

@epostma Thanks! Unfortunately, I think that it'll require that kernel work because according to the Maple Programming Guide, section 9.6 "Overloading Built-in Routines," select, etc., aren't on that short list of overloadables. This is surprising because I'd expect to be able to use select on any container object.