Double integral solver in TI-84
Here's the important part:

0→V
For(K,0,N-1
    A+KD+.5D→X
    V+fnInt(abs(Y₀-Y₉),X,A+KD,A+(K+1)D)fnInt(Y₈,Y,Y₀(X),Y₉(X))/(Y₀(X)-Y₉(X→V
End


I suggest you replace the fourth line with:

V+DfnInt(Y₈,Y,Y₀,Y₉)→V


Here's the reason: You should have thin enough slices that the volume in a single slice of width D can be approximated by DfnInt(Y₈,Y,Y₀,Y₉); that's more important than getting the area exactly right. This change will increase your overall speed by about 3x.

Once you've done that, you can drop the For( loop entirely, and use sum(seq( instead. This only helps the performance a little; your bottleneck is still the fnInt(. However it's better style.

Another trick: since your precision is limited more by the Midpoint rule calculation along the x axis than the accurate fnInt( along the Y axis, you may want to increase the tolerance of fnInt( to, say, .005. Then your entire loop becomes:

Dsum(seq(fnInt(Y₈,Y,Y₀,Y₉,.005),X,A+.5D,B,D→V


(This uses the optional 'step' argument to seq(.

Another possibility is to use the Midpoint rule on both axes instead of fnInt(. I'm not sure how much this would help speed or decrease accuracy, but I suspect the accuracy wouldn't be too much worse than using the Midpoint rule on one axis.